How to approximate gaussian kernel for image blur

The continuous Gaussian, whatever its dimension (1D, 2D), is a very important function in signal and image processing. As most data is discrete, and filtering can be costly, it has been and still is, subject of quantities of optimization and quantification/quantization schemes. In one 1D, the three most direct for a finite-length filter are illustrated below:

In solid blue, the continuous Gaussian. For an $L=7$-length, I choose a Gaussian with parameter  $s=\sqrt{L/2}$, in consistence with the first approximation $(P)$ given by the binomial formula or the Pascal triangle, as detailed by @Olli Niemitalo.
Each filter coefficient (red crosses) is of the form $$h_P[k] =  \frac{\binom{L}{k}}{2^L}\,.$$
Their sum is naturally equal to unity, and the coefficients are dyadic-rationals, in the form of $a/2^b$. They are practical for integer operations, on integer-valued pixels. As you can see, they are slightly above the Gaussian around the center, and below on the tails. When $L\to \infty, the shape tends to a certain bell function shape, see From Pascal's Triangle to the Bell-shaped Curve.
A second approximation (green circles), often called the truncated Gaussian, consists in sampling ($S$) $x\mapsto G(x)=\frac{1}{s\sqrt{\pi}}e^{-\frac{x^2}{s^2}}$ at integer values of $x$: $[-3,-2,\ldots,0,\ldots,3]$. They need to be normalized to 1-unit sum:
$$h_S[k] =  \frac{G(k-(L-1)/2)}{\sum_{l=1}^L G(l-(L-1)/2)}\,.$$
The third approximation ($A$), blue circles and bars) is area dependent: once again after normalization, the coefficients are proportional to the area under the Gaussian, on the interval $[k-1/2,k+1/2]$:
$$h_A[k] \propto \int_{k-1/2}^{k+1/2} G(x)dx\,.$$
The exact value depends on the quadrature formulae.
In our settings (with the chosen $s$), coefficients are, respectively from left to right (pascal, Sampling, Area):
$$
 \begin{Bmatrix}
0.0156 & 0.0232 & 0.0255\\
    0.0938 & 0.0968 & 0.0998\\
  0.2344 & 0.2282 & 0.2262\\
 0.3125 & 0.3036 & 0.2970\\
 0.2344 & 0.2282 & 0.2262\\
 0.0938 & 0.0968 & 0.0998\\
 0.0156 & 0.0232 & 0.0255\\
 \end{Bmatrix}
$$
Only the first column admits dyadic/integer operations. To work with integer arithmetic, one can round to the nearest integer. As $2^7=128$, this quantization choice (another approximation)  leads to:
$$\frac{1}{128}
 \begin{Bmatrix}
2 &  3 &  3\\
12 & 12 & 13\\
30 & 29 & 29\\
40 & 39 & 38\\
30 & 29 & 29\\
12 & 12 & 13\\
2 &  3 &  3\\
 \end{Bmatrix}
$$
With  $2^6=64$, you get:
$$\frac{1}{64}
 \begin{Bmatrix}
1 &  1 &  2\\
6 &  6 &  6\\
15 & 15 & 14\\
20 & 19 & 19\\
15 & 15 & 14\\
6 &  6 &  6\\
1 &  1 &  2\\
 \end{Bmatrix}
$$
In the case you mention in comments, the formula for $\sigma$ is a bit different: $\sigma=(3L+1)/20$, equal to $1.4$ for $L=7$. I don't know its origin yet. As said above, those are the main  principles. In the Gaussian Kernel Calculator demo, with those values, you get:
$$\begin{Bmatrix}
0.031251\\  0.106235    \\0.221252\\    0.282524    \\0.221252\\    0.106235\\  0.031251
\end{Bmatrix}
$$
which is quite close (up to the integral quadrature precision) to your:
$$\begin{Bmatrix}
0.03125\\0.109375\\0.21875\\0.28125\\0.21875\\0.109375\\0.03125
\end{Bmatrix}
$$
as the former, multiplied by 64, gives:
$$\begin{Bmatrix}
2.0001  \\   6.7990 \\  14.1601 \\  18.0815 \\  14.1601 \\   6.7990 \\   2.0001
\end{Bmatrix}
$$
The ($A$) approximation is possibly more accurate with respect to the  integral version of the convolution, with a
continuous Gaussian and a discrete or piecewise constant signal.
In 2D, one can tensorize two 1D filters, in which case the 2D filter has rank one. Or one
extends the above method to 2D.  The 2D coefficients can be obtained by two-dimensional integral estimations:
$$h_A[m,n] \propto \iint_{[m-1/2,m+1/2]\times  [n-1/2,n+1/2]} G(x,y)dxdy\,.$$
See for instance Gaussian smoothing, esp. Figure 3, Discrete approximation to Gaussian function with $\sigma=1.0$:


Effect of noise on high and low frequency components

This is actually due to the characteristics of the actual signal, rather than noise.

You might have heard of pink noise and its spectral characteristics $$S(f)\propto \frac{1}{f^\alpha}$$
Now forget about the term "noise" in here, and imagine a stochastic process whose spectral density vanishes following a $\frac{1}{f^\alpha}$ law, in which $\alpha$ is a constant.

Interestingly, this is the law that nature follows a lot. More specifically, many natural signals follow this law. It has been observed and actually been verified to be valid for instance in human voice, natural images, music, quasar light emissions, heart beat, firings of neurons, and even in financial markets. Note that natural is a key term here and by that, we refer to attributes of a natural environment. For a more detailed discussion see here and here. For 2D (image) signals we should expect $\alpha=2$.

Since according to the $\frac{1}{f}$ law the contribution of the actual [natural] signal will diminish at higher frequencies, we would expect more contribution of noise at these frequencies. Consequently, noise is usually dominant at high frequencies and the signal is more dominant at lower frequencies.


Where did the k of $a_k$ disappear from Fourier Reverse Transform if $\omega=\omega_0k$?

It is because $\omega_o$ approaches 0 such that $k\omega_o$ becomes a continuous function.

So the discrete function $X(jk\omega_o)$ becomes the continuous function $X(j\omega)$.

This is intuitive if you consider the Fourier Series expansion in that each component of the spectrum is spaced in frequency by $1/T$ Hz, (or $\omega_o = 2\pi/T$ radians/sec). Each of these components are given as $k\omega_o$ with $k$ being the set of all integers. So as $T$ gets larger, $1/T$ gets smaller, and hence $\omega_o$ gets smaller and as $T$ approaches infinity, $k \omega_o$ approaches a continuous function.

The concept behind the Fourier Series Expansion is that any continuous time function defined over time duration $T$ can be decomposed into a sum of sinusoids that each have a fundamental period of $T$ and all higher harmonics such that there are an integer number of cycles in $T$, the the waveform can repeat for all time and start and end the same way at all boundaries $t= NT$ for all integers $N$. For this reason no other frequency components can exist and the non-zero frequency spectrum is therefore discrete.

The Fourier Series Expansion given above in terms of sines and cosines is much more compact using exponential frequency notation (and ultimately much more intuitive once you understand that $e^{j\Omega t}$ is a phasor rotating counter-clockwise on the complex plane at rate $\Omega$ radians per second) :

$$x(t) = \sum_{n=-\infty}^{\infty}c_n e^{jn\Omega_o t}$$
Where

$e^{jn\Omega_o t} = \cos(n\Omega_o t) + j\sin(n\Omega_o t)$

$c_n = a_n - jb_n$

And $\Omega_o$ was used in place of $\omega_o$ just to match the plots I already had (where I was distinguishing radian frequency between analog and digital domains).

This becomes the Fourier Transform when we all $T$ to extend to infinity: the spacing given by $1/T$ between each of the components goes to zero, and the frequency spectrum therefore becomes continuous.

An easy way to observe this effect is by simply zero-padding a time domain waveform prior to taking the FFT: as you add more zeros, meaning $T$ gets larger, more and more samples will be interpolated in frequency as the spectrum approaches a continuous function.


Understanding a simple case of fft2 in python

When you construct re_field, ditch the np.abs. It turns all negative values of the sine wave positive. Instead, do a np.realon the result of the ifft2 to get rid of the imaginary parts caused by rounding errors. Then you will get a result perfectly matching the original.

In the magnitude spectrum you cannot see the sine wave you inserted, because its frequency is to low. The vertical black line is the DC part in x direction, which is zero, the horizontal white line is the DC part in y direction, which has full ampitude because there is only DC in y direction.

Try some higher frequencies, then you will get a better idea.


MATLAB: get members of histogram populations EXCLUDING those inbetween

Assumption

Your noise is Gaussian.

Proposed solution

Compute the histogram of your data. There will be two peaks, as in the histogram you've shown. Extract the approximate positions of the peaks and use these values as starting values to fit two Gaussian curves to both peaks, each with the fit parameters $\mu$ (mean signal, i.e. position of the Gaussian) and $\sigma$ (width of the Gaussian).

Your signal will then be the difference in positions of the Gaussians, and your noise is given by their widths.

Yes, there will still be contributions from the sides of the rect-pulse. However, there probably is no way to fully get rid of them, since you need some criterion to tell apart plateau values from values from the side lobes, what always becomes problematic if a point is close to the baseline/plateau. You could mitigate this effect by using the histogram values not only as datapoints for the fit, but also as fit weights. By doing so, the (low number of) side lobe values gets less weight compared to the (high number of) plateau/baseline values.


Question About Kailath's Paper - An Innovations Approach to Least Squares Estimation Part I: Linear Filtering in Additive White Noise

They used to have different convention for writing this stuff back then.
But actually what you saw is really simple.

It's all based on the Orthogonal Principle of MMSE.

They say in (9) that the Additive Noise is uncorrelated (By defining its Auto Correlation by Delat Function).
In (10) they say the optimal estimation of $ x \left( t \right) $ is given by $ \hat{x} \left( t \mid t \right) $ is a linear combination of $ y \left( t  \right) $ and $ z \left( t \right) $.
In (11) they exactly use the Orthogonal Principle which the estimator must obey.

In (12) they just derive the Correlation Matrix between the processes.
Since the linear estimator represent the correlation it is not surprising to see that it is the linear combination defined by the optimal filter $ g \left( \cdot \right) $.

The trick here is the fact (Using article notations) $ \overline{ \left( x \left( t \right) - \hat{x} \left( t \mid t \right) \right) v' \left( s \right) } = 0 $  which means (By Linearity of the Expectation) $ \overline{ x \left( t \right) v' \left( s \right) } = \overline{ \hat{x} \left( t \mid t \right) v' \left( s \right) } $.
Expand the right hand term by (10) by plugging $ \hat{x} \left( t \mid t \right) $ and you get (12).

By the way, Kaliath has much better written, in my opinion, paper on the subject called A View of Three Decades of Linear Filtering Theory.


I don't really understand what are you trying to achieve with vector B that has 6 elements. Pre-emphasis filter is defined as:

$$y[n]=x[n]-\alpha x[n-1]$$

Where $\alpha$ is $0.97$ in your case.

You can see that there are only two elements in vector B, right? So the code should be something like that:

B=[1, -0.97];
xd = filter(B, 1, data, [], 2);

You can verify the frequency response:

freqz(B, 1)


Duration of unknown rectangular pulse with additive white Gaussian noise

You want a method that removes noise while preserving edges. This cannot be achieved well by linear filtering, as you noticed yourself. I know of two approaches that might work well for your problem. The first is median filtering, where samples inside a window are replaced by their median. The following plot shows the result of median filtering with a window length of 25 samples (in red):

The other, more complex, approach is Total variation denoising, which works very well for piecewise constant signals. There is a very good description of total variation denoising including Matlab code available: link.


While I do not have proof, it seems pretty clear, given the "Characteristics of Compass signals reported as of May 2008 compared to GPS-L1CA" table at the Wikipedia page that you offered as an example, that the "(n)" notation refers to the chip rate.  You will notice that there is one BPSK(1) signal that has a chip rate of 1.023 Mchips/s.  There are four BPSK(2) signals listed, and all of them have a chip rate of 2.046 Mchips/s, which is twice 1.023 Mchips/s.  There are also three BPSK(10) signals that have a chip rate of 10.23 Mchips/s which is, of course, 10 times 1.023 Mchips/s.

It appears that the "(n)" notation means that the signal has a chip rate of $n*f_o$, and in the case of the Compass signals $f_o$ is 1.023 Mchips/s.


If the noise is extremely narrow-band (a perfect sinusoid), exactly integer periodic within the length n samples, and completely disjoint in the frequency domain from the non-noise signal of interest, then you might be able to zero FFT bin number (F_noise * n / F_sample_rate) and IFFT the result.

If the noise is wide band, or of a frequency not exactly at an FFT result bin center, then a time domain filter might be more suitable from removing noise without serious distortion artifacts.  Although you might be able to use an FFT/IFFT overlap add/save algorithm method to implement a time-domain FIR filter.


When input data into a low pass filter, are the first several data unusable?

You have to judge yourself if for your purposes you need to get rid of the first few output values. The phenomenon you observe is determined by two factors. The first is the delay of the filter (which is usually frequency dependent). This delay is a consequence of the causality of the filter. For a linear phase FIR filter that delay is independent of frequency and equals $(N-1)/2$ samples, where $N$ is the filter length (as pointed out in a comment by MBaz). The other factor lies in the nature of low pass filters and has nothing to do with a non-ideal (causal) implementation. It is the non-zero rise time that any low pass filter exhibits, even an ideal (non-causal) brick-wall filter. The rise time is proportional to the inverse of the cut-off frequency.

For minimum-phase filters with a very low order (as in your example), the first factor plays only a small role, so the delay you observe is mainly caused by the non-zero rise time of the low pass filter. If you decrease the cut-off frequency you will observe an increase in rise time and vice versa. Note that the input signal that your filter sees is actually a jump from zero (no input signal) to a DC value (plus some small signal), so you pretty much see the system's step response.

In general the response to a change in the input signal is called transient response, and it is not restricted to a change in DC value.


Play with a Gaussian Random Set in the Frequency Domain to Obtain Desired Effect in the Time Domain

You're asking to do a localized operation in time using the Frequency Domain.

It's going to be not elegant, really.

Here what you can do:

Define the input signal in frequency domain as $ X \left[ k \right] $.
Multiply by Linear Phase in Frequency Domain to apply circular shift in time domain.
This means samples 1-90 will become 11-100 and 91-100 will become 1-10. Lets' call this multiplied signal $ \hat{X} \left[ k \right] $.
Pre calculate the DFT of a signal which is 1 for the first 10 samples and then zero for samples 11-100. Lets call it $ H \left[ k \right] $. This is basically a windows on the first 10 samples. In order to extract the 10 first samples in time domain we would multiply by it. In the frequency domain it means we will apply Convolution. The result of the convolution will be a frequency domain signal which is equivalent of the DFT of only the 10 first samples padded with 90 zeros. Keep this result on side. Let's Call it $ Y \left[ k \right] = \hat{X} \left[ k \right] \ast H \left[ k \right] $.
Generate noise on frequency domain which is with the properties needed of generate what's needed on time domain. Convolve the generated noise with $ H \left[ k \right] $ which will again extract a signal which is equivalent to only 10 first samples in time domain. From this signal, lets call it $ Z \left[ k \right] $ subtract the signal $ Y \left[ k \right] $.
Add the result to $ \hat{X} \left[ k \right] $. Now the signal $ \hat{X} \left[ k \right] + Z \left[ k \right] - Y \left[ k \right] $ is what you're after.

Pay attention that doing this in Fourier domain might be much less efficient than doing it in time domain.
As you need to apply convolution twice instead of doing 2 DFT's.


Why are these results different about 'frequency shifting'?

These angular errors are in the order of $10^{-10}$. These are numerical inaccuracies caused by the fact that exp() can't be arbitrarily exact.

"Exact" would mean that exp(1i*2*pi*Fs*t) would just yield a vector [1 1 1 … 1]; if you'd hit that, the input to the fft would be identical, because multiplication with $1$ doesn't do anything.

So, always consider the error relative to the phenomenon you're observing. You wondered about something that was about 200 dB lower in amplitude than the observed quantity.


How to find spectral density of a signal whose correlation depends on time?

Your process is not stationary. As you already correctly noted, your autocorrelation function depends on $t$ and $\tau$. Let me call it $\varphi(\tau,t)$.

There are multiple ways of dealing with such cases. One is to simply consider Fourier transforms with respect to each of the time variables, treating them independently: The transform with respect to $\tau$ gives you frequency (say, $f$), where as the transform with respect to $t$ gives you a rate of change as in how fast do your statistics change, the latter often being referred to as Doppler frequency (say $\alpha$).

Now you can define four functions:

Time-varying ACF $\varphi(\tau,t)$
Time-varying Power spectral density $\Phi(f,t)$
Delay/Doppler cross spectral density $\varphi(\tau,\alpha)$: often also called scattering function
Frequency/Doppler power spectrum $\Phi(f,\alpha)$

These are also called the second set of Bello functions, the concrete naming of each of them varies widely across sources.

Another way of attacking the problem is to go to the Wigner-Ville distribution and its variants, have a look here.


Integer and a half periods of sine wave in DFT window

That's not true. Take a simple example:

N = 8; n = 0:N-1;
w1 = pi/4; w2 = 3*pi/8;
s1 = sin(n*w1); s2 = sin(n*w2);
S1 = fft(s1); S2 = fft(s2);
subplot(2,2,1),stem(n,s1)
subplot(2,2,2), stem(n,s2)
subplot(2,2,3), stem(n,abs(S1))
subplot(2,2,4), stem(n,abs(S2))


Optimally approximating the sign function by functions with compactly supported Fourier transform

This problem is in a way the dual of the problem of approximating a Hilbert transformer by a filter with finite memory. The frequency response of an ideal Hilbert transformer is
$$H(\omega)=-j\,\textrm{sgn}(\omega)\tag{1}$$
and the corresponding ideal impulse response is
$$h(t)=\frac{1}{\pi t}\tag{2}$$
We can approximate $(1)$ and $(2)$ by multiplying the ideal impulse response by a window of finite length:
$$\tilde{h}(t)=h(t)w(t)\tag{3}$$
where the window function $w(t)$ is a real-valued even and smooth function with finite support.
In your case you could use a window function $W(\omega)$ in the frequency domain to approximate a sign function:
$$f(t)=\textrm{p.v.}\frac{1}{j\pi}\int_{-\pi}^{\pi}\frac{W(\omega)}{\omega}e^{j\omega t}d\omega=\textrm{p.v.}\frac{1}{\pi}\int_{-\pi}^{\pi}\frac{W(\omega)}{\omega}\sin(\omega t)d\omega\tag{4}$$
Note that this is just a generalization of what you already did. You used a rectangle and a triangle for $W(\omega)$. However, there are many more "nicer" window functions that are smooth, and that may result in better approximations of the sign function.


How to Improve Least Means Squares (LMS) / NLMS Filter Performance?

If you want to implement the "standard" NLMS algorithm without cutting any corners, then you're probably not going to find a structure that is significantly more efficient. Block forms of LMS filtering aim to use fast convolution techniques (like overlap-save or overlap-add) to speed that part of the process. However, as you noted, the filter coefficients are only updated per block, as the filter must be constant over the block to use the fast convolution approach.

The highly recursive nature of NLMS is going to limit you if you want to keep the sample-by-sample update characteristic. While the filtering action is non-recursive, the filter coefficients at time instant N are a function of the coefficients at time instant N-1, which limits your ability to speed the process by using parallelism or block-oriented computation. As in most cases, there isn't a free lunch: if you want pure NLMS, then you're best off just implementing that.


What is the least JPG-compressible pattern? (camera shooting piece of cloth, scale/angle/lighting may vary)

Noise is quite good (hard to compress), but it becomes grey when looking from far, becoming easy to compress. A good pattern would be kind of fractal, looking similar at all scales.

Well, there is fractal noise.  I think Brownian noise is fractal, looking the same as you zoom into it.  Wikipedia talks about adding Perlin noise to itself at different scales to produce fractal noise, which is maybe identical, I'm not sure:

I don't think this would be hard to compress, though.  Noise is hard for lossless compression, but JPEG is lossy, so it's just going to throw away the detail instead of struggling with it.  I'm not sure if it's possible to make something "hard for JPEG to compress" since it will just ignore anything that's too hard to compress at that quality level.

Something with hard edges at any scale would probably be better, like the infinite checkerboard plane:

Also something with lots of colors.  Maybe look at actual fractals instead of fractal noise.  Maybe a Mondrian fractal?  :)


Flat vs frequency selective fading in MC-CDMA?

It would help if you mentioned where your references for MC-CDMA and your channel model come from.

The short answer is that it depends on the level of abstraction you want. For a flat channel, you would multiply by a complex-valued gain. For a frequency-selective channel, I would typically form the time domain signal and convolve with an impulse response (generated from some channel model that reflects your target environment). Perhaps if you provide us with a reference, I can provide additional assistance.


Getting a part of the invert of Si function

There are two problems, I think: fitting a parabola to the higher part, and moving the parabola to the next part of the sequence. This will address fitting the parabola.

I'd just use a least squares technique to find the $a$, $b$, and $c$ parameters.

This shows the part of the paper of interest:


determine background noise level via cross correlation

If you want to get an estimate of the noise, do the cross-correlation, and then multiply your transmitted sound by your cross-correlation peak and then subtract it from your received sounds.  Make sure, of course, that you use the cross-correlation to line up the transmitted sound with when it was received.

What remains after you subtract is the noise with some amount of (hopefully small) error.


Is there any way that we can perform speech recognition without using Fourier transforms?

All recognition tasks (doesn't even have to be speech recognition) are reductions of a very high-dimensional signal (your speech recording's dimension is the number of audio samples!) to a low-dimensional signal.

As such, it is generally advisable to transform the input signal through an easy operation to a representation where the dimensionality can be more easily reduced. For speech, the frequency domain is such a representation – physically, only a very limited set of discrete frequencies is sufficient to make up speech.

Hence, having a frequency-domain transform (typically, in speech recognition) early on in the process makes a lot of sense: after the transform, only very few coefficients contain information, and need to be evaluated further.

Hence, since speech models are based on the physics of speaking, and that is a time-frequency physical problem, it's very likely that an algorithm that actively tries to stay away from frequency domain will not perform good.

Nowadays, however, many classification problems are solved using Deep Neural Networks; these basically are just a network of nonlinear functions with a constant factor and offset for each of them. Training such a Neural Network to do something useful is just finding the right factors and offsets. There's heuristics that allow us to do that automatedly at large scale.

Such networks can of course be trained directly on discrete time-domain speech signals, and might have good recognition rates after sufficient training has happened.

Notice that training such Neural Networks is not a trivial task in itself; it's usually an intentionally or inherently stochastic procedure, and might or might not lead to a functioning classifier. There's a lot of experience and luck involved building a good Neural Network classifier instead of just "something that does something close to what I want".

However, the size and hence cost, and the effort spent on training, and the chance that learning doesn't lead to a working result, can be influenced by the way the signal is preprocessed. Think of it like this: Imagine you'd have to learn a language by listening to someone through a bad phone connection, compared to having a good phone at each end of the line that does some error correction and sound enhancements.

A typical preprocessing step here would be some sort of frequency-domain transform. Simply because the result of that more compactly describes the speech content than the time representation, it'll also make the Neural Network training's job of finding important and unimportant parts of the signal easier.

So, to conclude:

Yes, you can avoid the extremely simple Fourier transform (or one of the many other frequency-domain representations, many of which are based on the Fourier transform), it would probably make your life much harder than necessary.


Octave: Finding the mean fundamental frequency given an audio signal

Your code appears to be looking for zero crossings of the rectified signal, which won’t work reliably if your signal has strong harmonics or overtones of the fundamental frequency within your filter’s passband.  Strong overtones of low frequencies is very common in audio data.

As a likely better alternative, look up pitch detection/estimation algorithms, such as autocorrelation, ASDF and AMDF


What is the filter to get a heavy metal film ilustration?

There are two techniques that can achieve similar results, but if you caught that in a recent paper then you are probably after Style Transfer.

It is based on "teaching" a neural network how to combine specific elementary features in one context (e.g. Renaissance paintings) and then applying them in a different context (e.g. some random picture).

The other technique is "Colour Transfer" (or "image re-colouring") and that one attempts to copy just the colour tones and overall "feel" from one image to another but without changing the structure of its content. Here is a good overview and a particular method of applying it.

Hope this helps.


Phase response of $H(f)=e^{-j2{\pi}ft_0}$

Wouldn't it be much easier to just compare the left-hand side and the right-hand side of

$$e^{-j2\pi ft_0}=|H(f)|e^{j\phi(f)}\tag{1}$$

to see that $|H(f)|=1$ and $\phi(f)=-2\pi ft_0$?


Extracting direction of dominant sinusoids from an image in terms of angles

The following should work

Perform FFT

$$F(u,v) = \mathcal{F}[f](x,y)$$

Multiply magnitude by 2nd order circular harmonic and sum

$$c = \sum_u \sum_v |F(u,v)|  e^{i 2 \phi}$$

where $\phi = \arg (u + iv)$

Angle is half of argument

$$\theta = \arg(c) / 2 $$

Here is some matlab code:

I = rgb2gray(imread('eNb2pye.jpg'));
F = fft2(I);
[~,~,r,th] = FFTMesh(size(I));
c = sum(sum(abs(F) .* exp(1i.*2.*th)));
direction = angle(c)/2

which uses the following function to general the spectrum u,v and radial coordinates

function [ux,uy,r,th] = FFTMesh(sizeM)
rows = sizeM(1);
cols = sizeM(2);
[ux, uy] = meshgrid(([1:cols]-(fix(cols/2)+1))/(cols-mod(cols,2)), ...
    ([1:rows]-(fix(rows/2)+1))/(rows-mod(rows,2)));
ux = ifftshift(ux);   % Quadrant shift to put 0 frequency at the corners
uy = ifftshift(uy);
% Convert to polar coordinates
th = atan2(uy,ux);
r = sqrt(ux.^2 + uy.^2);
end


How to find the inflection point in a noisy curve?

I think Jazzmaniac's suggestion concerning Savitkzy-Golay filters makes a lot of sense. However, I do not suggest to simply smooth your data using these filters. What many people are not aware of is the fact that there are differentiating Savitzky-Golay filters, which smooth and compute the derivative at the same time. Such filters are implemented in Matlab: check this

You can use the smoothed derivative to detect your desired starting point. It will still require some heuristics and tweaking, but I do not see any other more straightforward way to do it.


Whats is the difference between FIR/IIR filters and Chebyshev/Butterworth filters

Chebychev and Butterworth filters can be transformed with the bilinear transform into IIR filters.  The amplitude and phase responses will be warped according to the bilinear transform.  This warp mainly makes for a higher cutoff rate than the corresponding analog filter (handwavingly, it's because the antialiasing filtering and sampling already does part of the work) which tends to be ok for most applications.

In contrast, using a bilinear transform on a Bessel filter is quite pointless since its salient design feature is a near constant phase delay in the pass band, and the bilinear transform warp messes this up.  So there is a separate digital equivalent to Bessel filters that isn't related via simple bilinear transform.

FIR filters, in contrast, don't have straightforward connections to analog filter designs.  Their design procedures are quite separate and only suitable for digital filters.


Telling apart a right traveling non-periodic signal from a left-traveling non-periodic signal, with a background periodic signal

If you can sample at any point in time, this shouldn't be that difficult.  The first step is to collect as series of frames in time.
In each frame you should be able to find a best fit sinusoid, or set of harmonic sinusoids, to your periodic signal.  Since this signal is assumed to have the same parameters across time, except the phase values which would need to be time adjusted, you can find an "average" value across all your frames for the periodic signal.
Subtract the best fit periodic signal from your frames.  This should leave the other two behind.  From your drawing, it appears that you have regions where one hump will b e zero, and others where the other hump is zero.  Those are the locations to read your humps at.
For more accuracy, once you have your hump signals estimated, go back and subtract them from your frames.  This should leave the periodic signal behind.  Once again, find the best fit sinusoid assembly.  This time your results should be more accurate.
Rinse and repeat.


they are distributed randomly, what I know that every column has 4 zeros and other values either 1 or -1
compression for that image can be done efficiently, Compared with the image which is full of 1 and -1 without any zeros

Since your spatial representation is already very sparse, it's likely that the DCT will reduce the compressability.
Simple calculation:
naive storage
That's a matrix of $7\cdot 8 =56$ entries, each of which takes one of 3 states, so 2 bits per pixel, so 112 bits in total.
non-zero position + sign storage
We can calculate the column entropy rather easily:

For the non-zero positions, you have 8 choose 4 different options, $\binom84=70$; so the info in that is $-\log_2\left(\frac1{70}\right)$; You can store that in 7 bit without thinking much. Just write a table of all possible non-zero positions, and number that.
If we consider the sign of the nonzero entries as bits, that's 4 bits of information.

So, the entropy in each column is $\log_2(70)+4\approx 10.13$ bit.
So, per column you need 11 bits, so with 7 columns in total 77 bits, and would, without any compression, have represented your whole image.
If we look at the 7 columns as whole, that's $70^7< 2^{43}$ possibilities of non-zero positions, i.e. we can reduce the overall storage needs to $43+7\cdot4=71$ bits. Again, no entropy coding necessary.
DCT first, entropy coding later
If you do a DCT on each column, you'll need all 8 entries; every DCT bin can take at least 5 values, so it needs 3 bit each, so 24 bit per column, so 148 bits in total.
Since the $70\cdot16$ (see above) different possible input columns are equally likely, so are the $70\cdot16$ possible DCTs; the information in seeing any one of these thus is $-\log_2\left(\frac{1}{70\cdot16}\right)\approx10.13$ bit. Since they don't have any correlation, entropy coding can, at the very best case compress your seven columns to 70.90 bits, ie. 71 bits in storage. But that case can never be achieved – that would work with an infinitely large image, not one with but 7 columns.
conclusion
for your 8×7 image:

intuitive storage format that just saves the non-zero positions + signs: 71 bits
unachievable asymptotic compression = entropy of the source: 70.9 bits

So, in this case, there's an "obvious" storage format with negligible redundancy. Thus, entropy coding (Huffmann etc) makes no sense. And, since your original image is already very sparse, any transformation, but especially the DCT, only makes compression harder.
So don't do the DCT in this synthetic example.


If I consider that your colored lines are multivariate signals (in Red, Green, Blue), that might be related through local offsets and scalings (stretch and squeeze), I would suggest the method in Correlation based dynamic time warping of multivariate time series, Expert Systems with Applications, 2012.

It combines Dynamic Time Warping (DTW), a widely used technique for comparison of time series data, and Principal Component Analysis (PCA)  to address the correlation between the variables.

yet, you have to balance, in your metric, the 'shape' part and the 'color' part of the data. The associated quantities do not live in the same space.


The feature vector is just the histogram of how many times a feature from each cluster appeared in the image.


Downsampling an image by an integer factor

Downsampling an image reduces the number of samples that can represent the signal. In terms of frequency domain, when a signal is downsampled, the high-frequency portion of the signal will be aliased with the low-frequency portion. When applied to image processing, the desired outcome is to preserve only the low-frequency portion. In order to do this, the original image needs to be preprocessed (alias-filtered) to remove the high-frequency portion so that aliasing will not occur.

The optimal digital filter to remove the high-frequency portion (with the sharpest cutoff) is sinc function. The reason is that the Sinc function's frequency domain representation is a nearly constant 1 over the entire low-frequency region, and nearly constant 0 over the entire high-frequency region.

$$\text{sinc}(x)=\frac{\sin(\pi x)}{\pi x}$$

The impulse response of the sinc filter is infinite. Lanczos filter is a modified sinc filter which attenuates the sinc coefficients and truncates them once the values drop to insignificance.

However, being optimal in frequency domain does not imply being optimal in human eyes. There are upsampling and downsampling methods that do not obey linear transformations but produce better results than linear ones.

With regard to the statement about $n \times n$, it is important to keep in mind that during image sampling, the choice of coordinates correspondence between the high-resolution signal and the low-resolution signal is not arbitrary, nor is it sufficient to align them to the same origin (0) on the real or discrete number line.

The minimum requirement in the coordinates correspondence is that

Upsampling an image containing arbitrary random values by an integer factor, then downsampling by the same integer factor, should result in the same image with minimal change numerically.
Upsampling/downsampling an image consisting of just one uniform value, followed by the opposite operation, should result in an image consisting of the same value uniformly, with minimal numerical deviations.
Repeatedly applying pairs of upsampling/downsampling should minimize the shift in image content as much as possible.


Representing a positive 0-limiting continuous time domain function, in terms of spectral samples, outside of a window

Finite $X_k$ samples will inevitably lead to a periodic approximation. The question, as is, is unfortunately and quite interestingly, unanswerable.

Indeed, under the standard formalisms, the reconstructed signal $\hat x(t)$ will be peridic:
$$
F{X(f)|||(f)}=\sum^{\inf}_{i=-\inf} X(i\Delta f)e^{i2\pi i\Delta f t}\Delta f = \hat x(t)
$$

If we truncate the signal $x(t)$ to be time-bounded, and using a similar argument than the Paley Wiener theorem to the direct Fourier Transform, we will obtain a exponentially converging spectrum on $X(f)$ when $f$ tends to infinity.

Based on the last, a logarithmic approximation could be sketched, in which the $f_i$ are being preferred when closest to zero, and sparse for increasing values:
$$
\hat x(t)=\sum^{\inf}_{f_i=-\inf} X(f_i)e^{i2\pi f_i t}\Delta f_i
$$

And under such aproximation, we could let the first terms of the sum to be again, approximated with known spectrum shapes (as $sinc(f)$ as you said). But that is not a train of frequency samples. This is a different work, and perhaps a different question?.


Total bandwidth in a network calculation

There are many variables involved, but you may proceed as follows. Assume the bit rate in a link is $R_b$. First of all, you choose a constellation; let's say you go with 16-QAM, which transmits four bits per symbol. Now you have a symbol rate $$R_1=R_b/4.$$ Next, you choose an encoder that gives you the BER you need. This encoder will have a certain code rate $R_c$. After encoding, you now have a symbol rate $$R_2=\frac{R_1}{R_c}.$$ Next, you need to take into account the overhead: you need framing and training data. Let's say the overhead is $x\%$. Then your final symbol rate is $$R_3=(1+x)R_2.$$ Now you can calculate the bandwidth you need: $$B=\gamma R_3,$$ where $\gamma$ is an "excess bandwidth" over the theoretical $B=R_3$. The excess bandwidth is typically between 1 and 2.

In summary: the symbol rate is given by the data rate, the constellation, the encoder, and the overhead. The bandwidth is determined by the symbol rate and the excess bandwidht factor.

Note that each of these variables affects the others. For example, if you choose QPSK instead of 16-QAM, you'll have better BER, but $R_1$ will be larger. However, now you need an encoder with smaller rate, so $R_2$ will be smaller. And so on. So, you need to find a combination of variables that you can implement within your constraints and that gives you the BER and data rate you need.


That really depends on your requirements and on your exact definition of bandwidth.

Telephone range (300 Hz-3.4 kHz) has okay intelligibility but the quality of the voice is fairly compromised. For example it's hard to tell the difference between "s" and "f" sounds.

80 Hz to 8 kHz gives you very good quality, although it's still possible for the human voice to make sounds outside of this range. Someone like James Earl Jones can easily go don to 50 Hz and specially trained singers can go a lot lower. Some fun examples at https://www.youtube.com/watch?v=lkdMl1ITw5M

If you want all of it, you should go for the range of hearing which is about 20 Hz to 20 kHz.


The problem with the moving average is that the average is not robust to the presence of the outliers - so you would need a very large window size to "dilute" the outliers.

Try a non-linear filter instead, like a median filter: Apply a median filter on your signal - you would need a window size of at least 300 samples. Compute the difference between the original signal and the median-filtered version. If the difference is above a threshold, replace the signal by the median-filtered version.

Here is some scilab code that attempts to implement this suggestion.  The results are plotted here; it seems to work nicely on the fabricated data.

function sm = smooth(x,len)
    sm = filter(ones(1,len)/len,1,x);
endfunction

N = 6000;
x = 0.1*rand(1,N,'normal');
y = cumsum(x);
y = smooth(y,100);

clf
subplot(211);
plot(y)

Njumps = 20;

jump_indices = round(rand(1,Njumps)*N);

jump_length = 0;

y2 = y;

for idx = jump_indices,
    y2(min(N,idx:(idx+jump_length))) = y(min(N,idx:(idx+jump_length))) + 1;
    //plot(min(N,idx:(idx+jump_length)),y2(min(N,idx:(idx+jump_length))),'r')
    jump_length = jump_length + 1;
end

plot(y2,'g');

filter_length = 10;
y3 = y;
for k = 1:N,
    y3(k) = median(y3(max(1,min(N,(k-filter_length/2):(k+filter_length/2)))));
end

plot(y3,'k')

subplot(212);
plot(y-y2);
plot(y-y3,'r');


How to identify the beginning and the end of a signal using MATLAB

There are many sources on Detection Theory. I think there are two relatively simple approaches to your problem. From the example graph I can clearly see a distinct signal. if the desired detected signal has the same trend/shape a good detection technique would be in the direction of matched filtering.

(1) Matched filtering is basically correlating the test signal with the measured signal. if the signals are very similar to each other, matched filtering would return a high value on which you could detect that with a threshold. A good starting point for signal detection would be this lecture slide.

(2) However, that might already be a bit too difficult since you mentioned to be new to signal processing. Another approach when I look at the signal is seeing that the noise is very low when the signal is not present. I'm not sure if the example signals have these distinct amplitudes when the signal is not present (yellow at ~15, purple at ~-730). But if this amplitude is know, we could simply check when the signal is different from this steady value. Or we could calculate how much the signal changes by calculating the gradient (or simply the difference between samples) for example. If the gradient is high, the signal would be increasing at a higher rate.

Maybe start by calculating the difference between samples, i.e.: $s[i]-s[i-1]$ for each sample. I highly encourage to just play with the signals when starting with signal processing. Try different functions on the samples.


Fast Fourier Transformation using parallel processing coding for image processing

Well, you should have a look at FFTW as their implementation allows you to use more than one thread to perform the FFT. One thing that should note with this library is that multi-threading is not necessarly beneficial for all problems sizes as stated here. For a software application that I'm working on, for example, using multiple threads slows down things (the biggest FFT being on 4096 points).


How exactly do I plot a constellation Diagram?

Typically a constellation diagram is not just the raw RF values because that would be extremely disorganized. You would end up with a bunch of points in a gaussian distribution about the origin (thermal noise) and a "ring" around it that represents your signal at any particular phase during the sampling.
What you typically want to plot is the baseband signal, the RF after it is filtered, mixed against the carrier wave, and coherent (phase-locked). This will indicate which constellation values are being received with some noise (hence why in Figure 7 they are scattered about the data positions). I have no idea what is going on in Figure 8, but it's not a constellation diagram. It looks like something of the author's own invention to better analyze the given data.


How to prepare and plot unequally spaced, irregular data to a contour plot or similar with MATLAB

I found a possible solution with the help of some collegue.

When gridis the matrix containing an x and a z and a value column (each having the same number of elements):

[X,Y] = ndgrid(linspace(min(grid(:,1)),max(grid(:,1)),150),linspace(min(grid(:,2)),max(grid(:,2)),150));
Z = griddata(grid(:,1),grid(:,2),grid(:,3),X,Y,'cubic');
contourf(X,Y,Z)
axis equal
scatter(grid(:,1),grid(:,2),10,grid(:,3),'wx');

So basically

a grid is generated with a relativly fine resolution (150x150 points) on basis of the min/max limits of the original data
the original data is interpolated onto this grid with a cubic interpolation method
contourf plot is generated and axis is set equal
the original data points are scattered onto the contourf plot to make sure that the interpretation of the data is somehow more reliable and one can see at once that the areas are interpolated.

I realized in between that the real difficulty is the obstacle area, as is gets wierdly interpolated although there should be no data, as none of the interpolation methods is used to such boundaries or optimized to those.


SNR for the magnitude of a complex signal

The magnitude $|y[n]|$ follows a Rice (aka Rician) PDF. You can find analytical expressions for the average, variance and SNR in Voltage Signal-to-Noise Ratio (SNR) Nonlinearity Resulting From Incoherent Summations and in The Squaring-Loss Paradox. These two references further investigate the SNR of the addition of independent magnitudes.


Lowpass filter impulse response in frequency domain

Discrete time domain comes with its magics and one of the funniest of them out there is the result of the multiplication of an impulse response $h[n]$ by $(-1)^n$, which is equivalent to an alternating version of the original response.

The key to its understanding comes from a very little effort of converting $(-1)^n$ into a complex exponential and then invoking a fundamental DTFT property, namely:

if $x[n]$ has the DTFT of $X(e^{j\omega})$ , then $x[n]e^{j\omega_0 n}$ has the DTFT of $X(e^{j(\omega-\omega_0)})$

Therefore, since $(-1)^n = e^{j\pi n}$, then DTFT of  $$h_2[n]=(-1)^n h_1[n]$$ is $$H_2(e^{j\omega}) = H_1(e^{j(\omega - \pi)})$$

Finally the interpretation of $H_1(e^{j(\omega - \pi)})$ is as follows, given that $H_1(e^{j\omega})$ is an ideal lowpass filter of cutoff frequency $\omega_c = \pi/2$ then $H_2(e^{j\omega})$ is a shifted by $\pi$ to the right version of $H_1(e^{j\omega})$ which will be an ideal highpass filter with a cutoff frequecy of $\omega_c = \pi /2$. And the best way to see this will come from drawing the magnitude reponses on the paper...


Given only the zeros and poles, the impulse response is determined except for a constant $K_i$, which you don't show.

From your calculations, it is clear that the energy of the impulse response for all filters is different.

If you normalize all filters so that all impulse responses have the same energy (i.e. $\sum_n |h_i[n]|^2$ is the same for the 4 filters), then the 4 curves for partial energies will end up at the same value, and it will be clear that $H_1$ has its energy distributed earlier in time than the others.


How to correctly compute the EEG Frequency Bands with Python?

I believe there is a much simpler way to do this with numpy.fft.rfft and numpy.fft.rfftfreq.

In the below example, I have two seconds of random data between 0.0 and 100.0 sampled at 512 Hz. I'm plotting the band amplitude below, but if you wanted power, it would just be the square of the values.

I think the reason you saw overly large values for Gamma, was that your gamma range is much larger than the others, and you're taking the sum of all the values in that range.  This is fine if you are going to compare one EEG set's gamma value to another (relative amp/power), but not if you're just looking at one set and want to compare gamma to say alpha.

If you just want to know what EEG state the person is in, you could use np.max instead of np.mean.  Then whichever band has the max will tell you about the person's state.

import numpy as np

fs = 512                                # Sampling rate (512 Hz)
data = np.random.uniform(0, 100, 1024)  # 2 sec of data b/w 0.0-100.0

# Get real amplitudes of FFT (only in postive frequencies)
fft_vals = np.absolute(np.fft.rfft(data))

# Get frequencies for amplitudes in Hz
fft_freq = np.fft.rfftfreq(len(data), 1.0/fs)

# Define EEG bands
eeg_bands = {'Delta': (0, 4),
             'Theta': (4, 8),
             'Alpha': (8, 12),
             'Beta': (12, 30),
             'Gamma': (30, 45)}

# Take the mean of the fft amplitude for each EEG band
eeg_band_fft = dict()
for band in eeg_bands:
    freq_ix = np.where((fft_freq >= eeg_bands[band][0]) &
                       (fft_freq <= eeg_bands[band][1]))[0]
    eeg_band_fft[band] = np.mean(fft_vals[freq_ix])

# Plot the data (using pandas here cause it's easy)
import pandas as pd
df = pd.DataFrame(columns=['band', 'val'])
df['band'] = eeg_bands.keys()
df['val'] = [eeg_band_fft[band] for band in eeg_bands]
ax = df.plot.bar(x='band', y='val', legend=False)
ax.set_xlabel("EEG band")
ax.set_ylabel("Mean band Amplitude")


Is it worth recording 4K resolution from a 1080p camera input?

I do not know how the choice between 1080p and 4K affects the recording quality in Ecamm live.
But it makes sense to convert a high-bitrate 1080p video to high-bitrate 4K video before uploading to YouTube. In my experience, the resolution increase will significantly improve YouTube video quality. This is because YouTube will use a significantly higher re-encoding bitrate at 4K compared to 1080p. Only a part of the re-encoding bitrate increase is consumed by an increase in the baseline bitrate demand when going from 1080p to 4K. Much of it will remain available for improving the video quality.
Here is as an indirect example a busy frame of a high-bitrate 1080p video converted to high-bitrate 4K video and uploaded to YouTube, and played back using 1080p and 4K settings:

The blocky artifacts in the 1080p stream are much more prominent.


Why do some radio stations still use amplitude instead of frequency modulation?

Because (as you note) all the receivers would have to be replaced.
Because international agreements allot that part of the spectrum to AM.
Because an FM transmitter with 5kHz modulation needs 25kHz channel spacing to avoid problems with adjacent channels.
Because an FM transmitter with 5kHz modulation only has an audio bandwidth of 3kHz - which is far less than AM (depending on the standard used 5khz, 10kHz or 15kHz.)


How to determine the filter type based on the impulse response?

If you know for sure that the filter type is one of the four given ones, then the exercise is quite straightforward. Look at the zeros of the corresponding transfer function

$$H(z)=\sum_{n=0}^{N-1}h[n]z^{-n}\tag{1}$$

where $N$ is the filter length. The zeros at DC ($z=1$) and at Nyquist ($z=-1$) are given by

$$H(1)=\sum_{n=0}^{N-1}h[n]\tag{2}$$

and

$$H(-1)=\sum_{n=0}^{N-1}h[n](-1)^n\tag{3}$$

It should be obvious that a (classic) low pass has a zero at Nyquist but no zero at DC. For a high pass filter the opposite is the case. A band pass filter has zeros at DC and at Nyquist, whereas a band stop filter has no zeros at DC or at Nyquist.

PS: What you read about the zero of the impulse response is either non-sense, or you misunderstood what they meant.


How to achieve a smoothed frequency phase response with no outliers

If your signals are not perfectly integer periodic in the FFT aperture (width), the phase measurements of those narrow band signals will toggle between FFT result bins.  This phase toggling will look like phase noise, but is actually a windowing artifact.  It is due to any window edge discontinuity from any non-exactly-integer-in-aperture periodicity in the signal (plus the fact that all the basis vectors of a DFT are perfectly circular in aperture).

To stop this phase toggling (apparent phase noise), perform an FFTShift before the FFT.  The fftshift will rotate the edge discontinuity away from sample 0, and reference the DFT phase estimate to the center of your data window, where it is likely continuous (instead of to a phase discontinuity between the last sample and first sample of the data).

Added:

There are (at least) 2 ways to perform a time-domain fftshift.  The first is to circularly rotate the data vector 50% after windowing but before performing the FFT.  The 2nd way is to complex conjugate every other result bin (odd ones usually) after the FFT.  Due to the DFT shift property, these two methods produce identical results.  Note that these are in the opposite order from doing a frequency domain fftshift.


just getting started in Signal Processing - easy question

The symbol $z^{-1}$ refers to a delay of one sample interval, no matter how long this interval actually is. It comes from the $\mathcal{Z}$-transform: if $X(z)$ is the $\mathcal{Z}$-transform of a data sequence $x[n]$, then $X(z)z^{-1}$ is the transform of the delayed sequence $x[n-1]$.

The transfer function (i.e. the quotient of the $\mathcal{Z}$-transforms of output and input sequences) is given by

$$H(z)=\frac14\left(1+z^{-1}+z^{-2}+z^{-3}\right)\tag{1}$$

which in the time (sample) domain means that the output sequence $y[n]$ and the input sequence $x[n]$ are related by

$$y[n]=\frac14\left(x[n]+x[n-1]+x[n-2]+x[n-3]\right)\tag{2}$$

This is of course the average of the current and the past three data samples. If, as in your example, all data samples have the value $8$, then the output is

$$y[n]=\frac14\left(8+8+8+8\right)=8$$


Recursive version of DFT as presented in Cooley-Tukey paper

The Cooley Tukey method allows general factorizations $N=N_1N_2$. The easiest way to see this is to use index mapping. Consider the length-$N$ DFT of a sequence $x[n]$:

$$X[k]=\sum_{n=0}^{N-1}x[n]W_N^{nk}\tag{1}$$

with $W_N=e^{-j2\pi/N}$. Let's assume that $N$ can be factored as $N=N_1N_2$. Now use the following index mapping

$$n=n_1N_2+n_2,\quad n_1=0,1,\ldots, N_1-1,\quad n_2=0,1,\ldots,N_2-1\\
k=k_1+N_1k_2,\quad k_1=0,1,\ldots, N_1-1,\quad k_2=0,1,\ldots,N_2-1\tag{2}$$

Plugging (2) into (1) gives

$$X[k_1+N_1k_2]=\sum_{n_2=0}^{N_2-1}\sum_{n_1=0}^{N_1-1}x[n_1N_2+n_2]W_N^{(n_1N_2+n_2)(k_1+k_2N_1)}\tag{3}$$

The complex exponentials in (3) can be simplified as follows:

$$W_N^{(n_1N_2+n_2)(k_1+k_2N_1)}=W_N^{N_2n_1k_1}W_N^{N_1N_2n_1k_2}W_N^{n_2k_1}W_N^{N_1n_2k_2}=\\=
W_{N_1}^{n_1k_1}\cdot 1\cdot W_{N}^{n_2k_1}W_{N_2}^{n_2k_2}\tag{4}$$

because $W_N^{N_1}=W_{N_2}$, $W_N^{N_2}=W_{N_1}$, and $W_N^{N_1N_2}=W_N^N=1$. Plugging (4) back into (3) and rearranging terms gives

$$X[k_1+N_1k_2]=\sum_{n_2=0}^{N_2-1}\left[W_N^{n_2k_1}\left(\sum_{n_1=0}^{N_1-1}x[n_1N_2+n_2]W_{N_1}^{n_1k_1}\right)\right]W_{N_2}^{n_2k_2}\tag{5}$$

Equation (5) can be interpreted as follows: the result can be computed by computing $N_2$ length-$N_1$ DFTs of subsampled versions of the input sequence $x[n]$ (this is the inner sum in (5)), then multiplying the result by twiddle factors $W_N^{n_2k_1}$, and finally computing $N_1$ length-$N_2$ DFTs (the outer sum in (5)). This procedure is most easily visualized using matrices:

read the data row-wise into a $N_1\times N_2$ matrix
compute DFTs over all $N_2$ columns (this can be done in-place)
multiply the result element-wise with twiddle factors
compute DFTs over all $N_1$ rows (in-place)
read out the data column-wise to get the final length-$N$ result

This little octave/Matlab script illustrates the procedure:

N = 1536;          % non-power of 2 FFT length used in LTE
N1 = 512; N2 = 3;  % 512 x 3 = 1536
x = randn(N,1);    % some input signal

% generate N1 x N2 matrix of input data, reading them row wise
X = reshape(x,[N2,N1]).';

% DFT over each column
X = fft(X);

% Multiplication with twiddle factors
X = X.*exp(-1i*2*pi/N*(0:N1-1)'*(0:N2-1));

% DFT over each row
X = (fft(X.')).';

% read out column wise
X = X(:);


Why is mutual information symmetric but conditional entropy isn't?

No,because while mutual information tells you how much uncertainty can a random variable $X$ remove from another random variable $Y$, the conditional entropy of $X$ given $Y$ tells you how much uncertainty remained in X after using the information (given by $I(X;Y)$) that $Y$ gave it.

$H(X|Y)=H(X)-I(X;Y)$

That's why the conditional entropy depends on the value of the entropy before the observation and the mutual information isn't, because it is only the difference ($\delta$) between two entropy states, before and after the observation.


Calculating the gain for adjusting sound level to X db

You don't need to convert the voltages to dB, just calculate the required gain (as a factor) and then convert this to dB if necessary, i.e.

G = Vout / Vin                 # required gain (factor)

G_dB = 20 * log10(G)           # required gain (dB)

where

Vin = input voltage
Vout = output voltage

If you really want to work with dB throughout then you can do it like this (should give the same answer):

Vin_dB = 20 * log10(Vin)       # input voltage (dB)

Vout_dB = 20 * log10(Vout)     # output voltage (dB)

G_dB = Vout_dB - Vin_dB        # required gain (dB)

G = 10 ^ (G_dB / 20)           # required gain (factor)


MATLAB: Why my Pixel Indices are so different from that in Spatial Coordinates

It is about Expressing Image Locations.

There are two kinds of image location expression,

One and the most common one is the

pixel indices:

in row and column

The other one is Intrinsic Coordinate System:

express in x and y:


Calculating maximum signal frequency knowing its sampling rate?

The maximum possible frequency in the sampled signal is half the sampling rate, assuming that you are not bandpass sampling and/or there is an anti-aliasing filter. The question is whether the signal has content up to that rate or whether it stops at a lower frequency.

One way to determine this would be to take your sampled signal $x[n]$ and transform it to the frequency domain via FFT to get $X[k]$.

Then find the largest value of $k$ for which
$$
\left |X[k] \right| > t
$$
for $0 \le k \le \frac{K}{2}$ where $K$ is the length of $X$ and $t$ is some "threshold" below which you decide the signal level is effectively zero.


Why does increasing FFT length (narrowing bandwidth) not decrease Noise Power per bin?

A spectrum analyzer typically has more than one mode. One mode normalizes the frequency output such that a pure tone of amplitude “one” with frequency centered on a bin center has an amplitude of “one” in the spectrum display.

For a pure tone, nonrandom, centered on a bin frequency of amplitude “one” , the corresponding uniformly weighted  N sample DFT bin magnitude will be “N”, so all the bins are reduced by dividing by N.  (before magnitude squaring)

An actual Spectrum Analyzer will include a fudge factor correction for nonuniform windows, and average multiple DFTs .

I haven’t seen too many spectrum analyzers that don’t have a DFT under the hood.  The last analog swept frequency analyzer I used was a very long time ago. I have seen a few modern versions at frequencies where converters aren’t good enough but DFT versions dominate most of the market.


If a square wave is a sum of odd harmonic impulses, why is it continuous in the frequency domain?

As @Hilmar mentioned I think you get confused between Square wave and  Rectangular function.
In Wikipedia about Square Wave :

A square wave is a non-sinusoidal periodic waveform in which the
amplitude alternates at a steady frequency between fixed minimum and
maximum values, with the same duration at minimum and maximum.

Which its Fourier Transform is only at harmonic frequencies and its value is equal to Fourier series coefficient.
but Rectangular Function:

The rectangular functionis defined as:

which look like this:

and its Fourier Transform is a sinc function that looks like this:


Obtain transfer function from signal flow graph

I am afraid that for complex systems finding such free standalone tool would be hard. Obviously since you mentioned MATLAB, you should use Simulink then or this tool:

Fully Automated Flowgraph Analysis Tool. Manual to be found here.

Alternatively you might look into this Mason.m tool. Or this masonry tool.

I know that these guys also did similar work once.


What transform can I use to "zoom" in a specific portion of a signal's spectrum?

You're probably referring to the zoom FFT. It's essentially a technique that allows for complexity reduction in the case where you have a small portion of a larger band that you'd like to analyze at high spectral resolution. It prevents you from having to calculate the high-resolution frequency content in the bands that you don't care about. Roughly, the algorithm can be summarized as follows:

Apply a bandpass filter around the region of interest, eliminating the components outside the band that you care about.
Decimate the signal by a factor $D$, such that the resulting sample rate still meets the Nyquist criterion for the filter's passband width. Depending upon where the band's center frequency was, this process might also involve frequency-translating the signal to baseband.
Perform a DFT on the signal. In order to get the same frequency resolution at the output, the "zoomed" transform only requires a transform length that is $\frac{1}{D}$-th of what you would need to use if you used the original, unfiltered signal for your analysis.


Where does this star came from? [EDIT - Detailed Question]

This seems as a starred transformation. It is a discrete version of the Laplace transformation. Laplace transformation is continuous by nature. Starred transformation is discrete, and the star denotes difference between continuous and analog domain.

This is commonly used when teaching about the Z transform. It tends to be confusing, so the teachers typically just skim over it.

Here are some explanation about it:  https://en.wikipedia.org/wiki/Starred_transform

www.fixoncloud.com/Home/LoginValidate/OneProblemComplete_Detailed.php?problemid=934


The formula you need, which I don't see in the scanned pages, is the N-dimensional Gaussian distribution of the received vector given that the signal $s_m$ was transmitted:

$$f(\vec{r}|s_m)=\frac{1}{(\pi N_0)^{N/2}}\exp\left[-\frac{\sum_{j=1}^N(r_j-x_{mj})^2}{N_0}\right]\tag{1}$$

The maximum likelihood receiver seeks to maximize this conditional probability by choosing the appropriate $s_m$. Since the logarithm is a monotonic function, you can as well maximize the logarithm of (1):

$$\log(f(\vec{r}|s_m))=-\frac{N}{2}\log(\pi N_0)-\frac{1}{N_0}\sum_{j=1}^N(r_j-x_{mj})^2=\\=-\frac{N}{2}\log(\pi N_0)-\frac{1}{N_0}||\vec{r}-s_m||^2\tag{2}$$

Since the first term in (2) is constant (i.e. independent of $s_m$), maximizing (2) is equivalent to minimizing the distance $||\vec{r}-s_m||^2$. If $||\vec{r}-s_m||^2$ is minimal, then obviously also $||\vec{r}-s_m||$ is minimal. So we arrive at the ML decision rule given in your book, i.e. choose the signal $s_m$ for which the distance to the received signal is a minimum.


Audio Activity Detection and Classification

Here is an recent introduction video on the subject (2019): Audio Classification with Machine Learning


Examples of Independent and uncorrelated data in real-life, and ways to measure/detect them

A few elements... (I know that this is not exhaustive, a more complete answer should probably mention moments)

Q1

To check whether two distributions are independent, you need to measure how similar their joint distribution $p(x,y)$ is to the product of their marginal distribution $p(x) \times p(y)$. To this purpose, you can use any distance between distributions. If you use the Kullback-Leibler divergence to compare those distributions, you will consider the quantity:

$\int_x \int_y p(x, y) \log \frac{p(x, y)}{p(x) p(y)} dx dy$

And you will have recognized... the Mutual Information! The lower it is, the more independent the variables are.

More practically, to compute this quantity from your observations, you can either estimate the densities $p(x)$, $p(y)$, $p(x, y)$ from your data using a Kernel density estimator and do a numerical integration on a fine grid ; or just quantify your data into $N$ bins and use the expression of the Mutual Information for discrete distributions.

Q2

From the Wikipedia page on statistical independence and correlation:

At the exception of the last example, these 2D distributions $p(x, y)$ have uncorrelated (diagonal covariance matrix), but not independent, marginal distributions $p(x)$ and $p(y)$.

Q3

There are indeed situations in which you might look at all the values of the cross-correlation functions. They arise, for example, in audio signal processing. Consider two microphones capturing the same source, but distant from a few meters. The cross-correlation of the two signals will have a strong-peak at the lag corresponding to the distance between microphones divided by the speed of sound. If you just look at the cross-correlation at lag 0, you won't see that one signal is a time-shifted version of the other one!


What is the difference between t=n*Ts sampling vs impulse train sampling?

I cannot recall a meaningful difference between the ideal impulse train sampling and the uniform sampling relation indicated by the expression $t_n = n T_s$.

Given a continuous-time bandlimited signal $x_c(t)$, when you sample this with an ideal impulse train $\delta_{T_s}(t) = \sum_k \delta(t - k T_s)$, then the relation between the obtained discrete-time sequence $x[n]$ and the continuous-time signal will be given by :

$$ x[n] = x_c(t_n) = x_c(n T_s)$$

Hence, the latter expression is a consequence of the former operation.


Why do cutoff frequencies of a high-pass digital filter differ from the ones i used to transform it from a low-pass analog filter in MATLAB?

The bilinear transform warps the frequency axis, that's why cut-off frequencies of the analog filter are not necessarily mapped to the correct frequencies in the digital domain. You can use the prewarping option by providing an optional parameter to bilinear.m. Check the corresponding mathworks documentation.


Your answer is right, assuming you posted the question right.

But you better use the standard notation as Dilip Sarwate already indicated; $u[n]$ is the unit-step and $\delta[n]$ is the unit impulse. Then

$$ 0.5^n \delta[n-2] \star \delta[n] = 0.5^2 \delta[n-2] = \begin{cases} { 0.25 ~~~, ~~~n= 2 \\ 0.00 ~~~,~~~n \neq 2 } \end{cases}
$$

you can get the answer for the second case, exactly in the same way.


The first term is just a constant, so it's not relevant to periodicity. The sequence as written is periodic with a fundamental frequency of 1/6.

Your cited solution makes no sense to me.


You are correct that $y[n]$ only depends on past input values for $n \gt 0$, but when $n \lt 0$, that is not the case --- then $y[n]$ is dependent on future input values, regardless of what $x$ is.

To see this:
$$
y[-2] = x[-(-2)] = x[2]
$$
so the output at time $n=-2$ is dependent on the input at time $n=2$.


What Approximation Techniques Exist for Computing the Square Root?

if you want a cheap and dirty optimized power-series expansion (the coefficients for Taylor series converge slowly) for sqrt() and a bunch of other trancendentals, i have some code from long ago.  i used to sell this code, but no one has paid me for it for nearly a decade.  so i think i'll release it for public consumption.  this particular file was for an application where the processor had floating point (IEEE-754 single precision) and they had a C compiler and dev system, but they did not have (or they didn't want to link in) the stdlib that would have had the standard math functions.  they did not need perfect precision, but they wanted things to be fast.  you can pretty easily reverse engineer the code to see what the power series coefficients are and write your own code.  this code assumes IEEE-754 and masked off the bits for mantissa and exponent.

it appears that the "code markup" that SE has is unfriendly with angle characters (you know ">" or "<"), so you'll probably have to hit "edit" to see all of it.

//
//    FILE: __functions.h
//
//    fast and approximate transcendental functions
//
//    copyright (c) 2004  Robert Bristow-Johnson
//
//    rbj@audioimagination.com
//

#ifndef __FUNCTIONS_H
#define __FUNCTIONS_H

#define TINY 1.0e-8
#define HUGE 1.0e8

#define PI              (3.1415926535897932384626433832795028841972)        /* pi */
#define ONE_OVER_PI     (0.3183098861837906661338147750939)
#define TWOPI           (6.2831853071795864769252867665590057683943)        /* 2*pi */
#define ONE_OVER_TWOPI  (0.15915494309189535682609381638)
#define PI_2            (1.5707963267948966192313216916397514420986)        /* pi/2 */
#define TWO_OVER_PI     (0.636619772367581332267629550188)
#define LN2             (0.6931471805599453094172321214581765680755)        /* ln(2) */
#define ONE_OVER_LN2    (1.44269504088896333066907387547)
#define LN10            (2.3025850929940456840179914546843642076011)        /* ln(10) */
#define ONE_OVER_LN10   (0.43429448190325177635683940025)
#define ROOT2           (1.4142135623730950488016887242096980785697)        /* sqrt(2) */
#define ONE_OVER_ROOT2  (0.707106781186547438494264988549)

#define DB_LOG2_ENERGY          (3.01029995663981154631945610163)           /* dB = DB_LOG2_ENERGY*__log2(energy) */
#define DB_LOG2_AMPL            (6.02059991327962309263891220326)           /* dB = DB_LOG2_AMPL*__log2(amplitude) */
#define ONE_OVER_DB_LOG2_AMPL   (0.16609640474436811218256075335)           /* amplitude = __exp2(ONE_OVER_DB_LOG2_AMPL*dB) */

#define LONG_OFFSET     4096L
#define FLOAT_OFFSET    4096.0

float   __sqrt(float x);

float   __log2(float x);
float   __exp2(float x);

float   __log(float x);
float   __exp(float x);

float   __pow(float x, float y);

float   __sin_pi(float x);
float   __cos_pi(float x);

float   __sin(float x);
float   __cos(float x);
float   __tan(float x);

float   __atan(float x);
float   __asin(float x);
float   __acos(float x);

float   __arg(float Imag, float Real);

float   __poly(float *a, int order, float x);
float   __map(float *f, float scaler, float x);
float   __discreteMap(float *f, float scaler, float x);

unsigned long __random();

#endif

//
//    FILE: __functions.c
//
//    fast and approximate transcendental functions
//
//    copyright (c) 2004  Robert Bristow-Johnson
//
//    rbj@audioimagination.com
//

#define STD_MATH_LIB 0

#include "__functions.h"

#if STD_MATH_LIB
#include "math.h"   // angle brackets don't work with SE markup
#endif

float   __sqrt(register float x)
    {
#if STD_MATH_LIB
    return (float) sqrt((double)x);
#else
    if (x > 5.877471754e-39)
        {
        register float accumulator, xPower;
        register long intPart;
        register union {float f; long i;} xBits;

        xBits.f = x;

        intPart = ((xBits.i)>>23);                  /* get biased exponent */
        intPart -= 127;                             /* unbias it */

        x = (float)(xBits.i & 0x007FFFFF);          /* mask off exponent leaving 0x800000*(mantissa - 1) */
        x *= 1.192092895507812e-07;                 /* divide by 0x800000 */

        accumulator =  1.0 + 0.49959804148061*x;
        xPower = x*x;
        accumulator += -0.12047308243453*xPower;
        xPower *= x;
        accumulator += 0.04585425015501*xPower;
        xPower *= x;
        accumulator += -0.01076564682800*xPower;

        if (intPart & 0x00000001)
            {
            accumulator *= ROOT2;                   /* an odd input exponent means an extra sqrt(2) in the output */
            }

        xBits.i = intPart >> 1;                     /* divide exponent by 2, lose LSB */
        xBits.i += 127;                             /* rebias exponent */
        xBits.i <<= 23;                             /* move biased exponent into exponent bits */

        return accumulator * xBits.f;
        }
     else
        {
        return 0.0;
        }
#endif
    }

float   __log2(register float x)
    {
#if STD_MATH_LIB
    return (float) (ONE_OVER_LN2*log((double)x));
#else
    if (x > 5.877471754e-39)
        {
        register float accumulator, xPower;
        register long intPart;

        register union {float f; long i;} xBits;

        xBits.f = x;

        intPart = ((xBits.i)>>23);                  /* get biased exponent */
        intPart -= 127;                             /* unbias it */

        x = (float)(xBits.i & 0x007FFFFF);          /* mask off exponent leaving 0x800000*(mantissa - 1) */
        x *= 1.192092895507812e-07;                 /* divide by 0x800000 */

        accumulator = 1.44254494359510*x;
        xPower = x*x;
        accumulator += -0.71814525675041*xPower;
        xPower *= x;
        accumulator += 0.45754919692582*xPower;
        xPower *= x;
        accumulator += -0.27790534462866*xPower;
        xPower *= x;
        accumulator += 0.12179791068782*xPower;
        xPower *= x;
        accumulator += -0.02584144982967*xPower;

        return accumulator + (float)intPart;
        }
     else
        {
        return -HUGE;
        }
#endif
    }

float   __exp2(register float x)
    {
#if STD_MATH_LIB
    return (float) exp(LN2*(double)x);
#else
    if (x >= -127.0)
        {
        register float accumulator, xPower;
        register union {float f; long i;} xBits;

        xBits.i = (long)(x + FLOAT_OFFSET) - LONG_OFFSET;       /* integer part */
        x -= (float)(xBits.i);                                  /* fractional part */

        accumulator = 1.0 + 0.69303212081966*x;
        xPower = x*x;
        accumulator += 0.24137976293709*xPower;
        xPower *= x;
        accumulator += 0.05203236900844*xPower;
        xPower *= x;
        accumulator += 0.01355574723481*xPower;

        xBits.i += 127;                                         /* bias integer part */
        xBits.i <<= 23;                                         /* move biased int part into exponent bits */

        return accumulator * xBits.f;
        }
     else
        {
        return 0.0;
        }
#endif
    }

float   __log(register float x)
    {
#if STD_MATH_LIB
    return (float) log((double)x);
#else
    return LN2*__log2(x);
#endif
    }

float   __exp(register float x)
    {
#if STD_MATH_LIB
    return (float) exp((double)x);
#else
    return __exp2(ONE_OVER_LN2*x);
#endif
    }

float   __pow(float x, float y)
    {
#if STD_MATH_LIB
    return (float) pow((double)x, (double)y);
#else
    return __exp2(y*__log2(x));
#endif
    }

float   __sin_pi(register float x)
    {
#if STD_MATH_LIB
    return (float) sin(PI*(double)x);
#else
    register float accumulator, xPower, xSquared;

    register long evenIntPart = ((long)(0.5*x + 1024.5) - 1024)<<1;
    x -= (float)evenIntPart;

    xSquared = x*x;
    accumulator = 3.14159265358979*x;
    xPower = xSquared*x;
    accumulator += -5.16731953364340*xPower;
    xPower *= xSquared;
    accumulator += 2.54620566822659*xPower;
    xPower *= xSquared;
    accumulator += -0.586027023087261*xPower;
    xPower *= xSquared;
    accumulator += 0.06554823491427*xPower;

    return accumulator;
#endif
    }

float   __cos_pi(register float x)
    {
#if STD_MATH_LIB
    return (float) cos(PI*(double)x);
#else
    register float accumulator, xPower, xSquared;

    register long evenIntPart = ((long)(0.5*x + 1024.5) - 1024)<<1;
    x -= (float)evenIntPart;

    xSquared = x*x;
    accumulator = 1.57079632679490*x;                       /* series for sin(PI/2*x) */
    xPower = xSquared*x;
    accumulator += -0.64596406188166*xPower;
    xPower *= xSquared;
    accumulator += 0.07969158490912*xPower;
    xPower *= xSquared;
    accumulator += -0.00467687997706*xPower;
    xPower *= xSquared;
    accumulator += 0.00015303015470*xPower;

    return 1.0 - 2.0*accumulator*accumulator;               /* cos(w) = 1 - 2*(sin(w/2))^2 */
#endif
    }

float   __sin(register float x)
    {
#if STD_MATH_LIB
    return (float) sin((double)x);
#else
    x *= ONE_OVER_PI;
    return __sin_pi(x);
#endif
    }

float   __cos(register float x)
    {
#if STD_MATH_LIB
    return (float) cos((double)x);
#else
    x *= ONE_OVER_PI;
    return __cos_pi(x);
#endif
    }

float   __tan(register float x)
    {
#if STD_MATH_LIB
    return (float) tan((double)x);
#else
    x *= ONE_OVER_PI;
    return __sin_pi(x)/__cos_pi(x);
#endif
    }

float   __atan(register float x)
    {
#if STD_MATH_LIB
    return (float) atan((double)x);
#else
    register float accumulator, xPower, xSquared, offset;

    offset = 0.0;

    if (x < -1.0)
        {
        offset = -PI_2;
        x = -1.0/x;
        }
     else if (x > 1.0)
        {
        offset = PI_2;
        x = -1.0/x;
        }
    xSquared = x*x;
    accumulator = 1.0;
    xPower = xSquared;
    accumulator += 0.33288950512027*xPower;
    xPower *= xSquared;
    accumulator += -0.08467922817644*xPower;
    xPower *= xSquared;
    accumulator += 0.03252232640125*xPower;
    xPower *= xSquared;
    accumulator += -0.00749305860992*xPower;

    return offset + x/accumulator;
#endif
    }

float   __asin(register float x)
    {
#if STD_MATH_LIB
    return (float) asin((double)x);
#else
    return __atan(x/__sqrt(1.0 - x*x));
#endif
    }

float   __acos(register float x)
    {
#if STD_MATH_LIB
    return (float) acos((double)x);
#else
    return __atan(__sqrt(1.0 - x*x)/x);
#endif
    }

float   __arg(float Imag, float Real)
    {
#if STD_MATH_LIB
    return (float) atan2((double)Imag, (double)Real);
#else
    register float accumulator, xPower, xSquared, offset, x;

    if (Imag > 0.0)
        {
        if (Imag <= -Real)
            {
            offset = PI;
            x = Imag/Real;
            }
         else if (Imag > Real)
            {
            offset = PI_2;
            x = -Real/Imag;
            }
         else
            {
            offset = 0.0;
            x = Imag/Real;
            }
        }
     else
        {
        if (Imag >= Real)
            {
            offset = -PI;
            x = Imag/Real;
            }
         else if (Imag < -Real)
            {
            offset = -PI_2;
            x = -Real/Imag;
            }
         else
            {
            offset = 0.0;
            x = Imag/Real;
            }
        }

    xSquared = x*x;
    accumulator = 1.0;
    xPower = xSquared;
    accumulator += 0.33288950512027*xPower;
    xPower *= xSquared;
    accumulator += -0.08467922817644*xPower;
    xPower *= xSquared;
    accumulator += 0.03252232640125*xPower;
    xPower *= xSquared;
    accumulator += -0.00749305860992*xPower;

    return offset + x/accumulator;
#endif
    }

float   __poly(float *a, int order, float x)
    {
    register float accumulator = 0.0, xPower;
    register int n;

    accumulator = a[0];
    xPower = x;
    for (n=1; n<=order; n++)
        {
        accumulator += a[n]*xPower;
        xPower *= x;
        }

    return accumulator;
    }

float   __map(float *f, float scaler, float x)
    {
    register long i;

    x *= scaler;

    i = (long)(x + FLOAT_OFFSET) - LONG_OFFSET;         /* round down without floor() */

    return f[i] + (f[i+1] - f[i])*(x - (float)i);       /* linear interpolate between points */
    }

float   __discreteMap(float *f, float scaler, float x)
    {
    register long i;

    x *= scaler;

    i = (long)(x + (FLOAT_OFFSET+0.5)) - LONG_OFFSET;   /* round to nearest */

    return f[i];
    }

unsigned long __random()
    {
    static unsigned long seed0 = 0x5B7A2775, seed1 = 0x80C7169F;

    seed0 += seed1;
    seed1 += seed0;

    return seed1;
    }


Looking for pratical quantitative comparison metrics for scaled, delayed and warped Signals

I'm answering the question the way I understood it - How can one find a similarity measure which isn't sensitive to scaling and shifting.

An approach could be borrowed from the Computer Vision world by comparing Shift and Scale Invariant features between the two signals.

I'm not sure it will work for measuring the quality of recovering signals but it certainly will tell two signals are similar given the hierarchy between their features is similar and the features themselves are similar.


This concept is related to the topological skeletonization of an object in shape analysis. It produces a thin sert of lines equidistant to the boundaries, preserving  geometrical and topological properties of the object. Some use the term medial axis. Some points possess a specific interest: branch or junction points (where sereval lines meet, i.e. a point with more than two neighbors) and termination, terminating or end points (a point with less than two neighbors), located at extremities.


Technically this is not a MATLAB-esque forum, but I can explain the steps in more detail for you: Suppose your input signal is $x[n]$, and its DFT is $X(f)$. For real signals you may use the one-sided DFT, since the other half would be redundant when you look at its Power Spectral Density. (PSD).

Once you compute the DFT of your signal, the PSD is simply $|X(f)|^2$. That is, you need to take the absolute magnitude of your DFT result, squared.

You now need to normalize the PSD such that it can be viewed as a Probability Density Function, (PDF). Thus, a normalized PSD, (let us call it $PSD_n$) will simply be:

$$
PSD_n(f) = \frac{PSD(f)}{\sum_{f=\frac{-fs}{2}}^{f = \frac{f_s}{2}} PSD(f)}
$$

Finally, your spectral entropy will be:

$$
E = -\sum_{f=\frac{-fs}{2}}^{f = \frac{f_s}{2}} PSD_n(f) log_2[PSD_n(f)]
$$


Natural images can seem complicated. However, the set of humanly interpretable images is relatively limited, with respect to all possible images ($24^{x\textrm{millions of pixels}}$).
There is a belief that their structure can be better summarized. Indeed, lossless compaction can achieve 2-3 factor ratios, and lossy compression about 8-16  with little perceivable loss.

Most of the present image compression methods (that I  know of) are based on some kind of interpolation/extrapolation schemes. The rationale is to predict some samples, and to code the hopefully smaller prediction errors either on less bits (lossless), or with less precision (lossy). Those schemes can be applied either in the original space domain, or in a transformed space: either on fixed bases (Fourier, discrete cosine transforms, wavelets, etc.), on adaptive bases (PCA, NMF, KLT), learnt dictionaries, machine learning or artificial intelligence dabases.

Adaptivity sounds sound. One can thus adapt to image features. Globally on images? Why not, but the computational prize can be important: eigen-algorithms can be $O(N^3)$. One can reduce the size to smaller image patches. You have to compute eigenvectors for each patch, passed on its second-ordre statistics (covariance).
Then, less memory is required, and this can be parallelized easily. Adaptive bases are thought optimal in several cases to pack information. But they come with a price for compression. Since bases are different from patch to patch, one should code the eigenvectors and transmit them as well. For traditional rate/distortion rates, this is not favorable.

So is there a good-enough fixed-basis? Among the many choices, it turned out than one discrete cosine transform, the DCT-II, could do a good job, either theoretically or practically. This is celebrated in a 1974 paper by N. Ahmed et al.: Discrete cosine transform
It is better suited to small patches, and possesses fast algorithms. There are other options, adaptations, but let's keep it simple for now.

So, assuming that one is working on patches, and computing covariance, DCT should be close to KLT in energy, and somehow better as one does not have to code vectors.

Below are previous answers related to properties of the DCT

How to understand the de-correlation property of DCT? What does de-correlation mean?
Could a DCT be used for an audio magnitude spectrum rather than DFT?
Why do we use scaling factors in DCT? Are they important?
Number of cosine waves in DCT
How to get the DCT energy distribution image in MATLAB,or we may call it the DCT coefficients distribution image?


Why is my bode plot showing periodic behavior?

When you are discretizing you make transformation $s \to e^{s T}$. This transforms the left half plane into unit disc. But actually it only transforms the $(-j\omega_s/2, j\omega_s/2)$ band and repeats itself by circling for outside of the band, where $\omega_s = 2 \pi / T$. This is because of the periodicity of exponential function.

When you plot the frequency response, you plot $|G(j\omega)|=|\hat{G}(e^{j\omega T})|$ for all $\omega \in [0, \infty)$ where $\hat{G}$ is the discrete transfer function. Therefore the frequency response of any sampled system is periodic.


Is it possible to estimate variance of noise for a step answer signal?

Hi: In order to estimate the variance, you need to have an underlying model for your signal. So, suppose that the model is

$y_{t+1} = y_t + \epsilon_t$ $~\forall ~ t = 1,\ldots n $.

assuming that $E(\epsilon_t) = 0$ and $var(\epsilon_t) = \sigma^2$.

In this case, you would difference your data in order to get estimates
of  $\epsilon_{t}$ at each time $t$, and call the estimates
$\hat{\epsilon}_t, ~\forall i = 2, \ldots n$.

Then, then an unbiased estimate of $\sigma^2 , \hat{\sigma}^2$, $= \frac{1}{n-2} \sum_{i=2}^{n}  (\hat{\epsilon_i} - 0)^2$.

But that estimate is critically dependent on the assumed model for the response $y_{t}$ ( this model is sometimes used in finance for the log price at time t+1 and is referred to as a random walk model ).  A different model will lead to a different estimate of the noise because the assumptions about the noise will be quite different in each case. I hope this helps some.


Why are analytic FFT and computational FFT not matching?

The fourier transform is defined for both positive and negative frequencies.  It looks like your first plot is just the positive frequencies, whereas the second plot contains both the positive and negative frequencies.  It's  common in computational fft routines to have the positive frequencies occupy elements 1 to N/2 in the results vector, and then elements N/2 to N are the negative frequencies.  Also, computational results are often unscaled.  You'll need to divide the results by the length of the window. Here is a short article I wrote explaining this for matlab http://mechanicalvibration.com/Making_matlab_s_fft_functio.html  Your results for python should be similar.


Your problem is called an assignment problem. Check out the Hungarian algorithm. Instead of average BER cost you can minimize total cost as it gives the same optimal solution.


Understanding spectra in polyphase filter banks

Since you define the input to the filter $p_1[n]$ as $x_1[n]=x[2n-1]$ and the channel impulse response as $p_1[n] = p[2n+1]$ the spectrum at the output of that channel will be as follows:

First express $X_1(e^{j\omega})$, the DTFT of the input to the channel no 1.

$$ x[n] \longrightarrow \boxed{z^{-1}} \longrightarrow v[n]=x[n-1] \longrightarrow \boxed{\downarrow2} \longrightarrow w[n]=x[2n-1] $$

Now as can be seen, $V(e^{j\omega}) = X(e^{j\omega}) e^{-j\omega} $ and also $W(e^{j\omega}) = \frac{1}{2} \sum_{k=0}^{1} V( e^{j( \frac{\omega-2\pi k}{2} )  })$ When combined it yields;

$$ X_1(e^{j\omega}) = \frac{1}{2} \sum_{k=0}^{1} X( e^{j( \frac{\omega-2\pi k}{2} )  }) e^{-j( \frac{\omega-2\pi k}{2} )} $$

In this case

$$ X_1(e^{j\omega}) = \frac{1}{2} e^{-j\frac{\omega}{2}} \cdot(X(e^{j\frac{\omega}{2}}) - X(e^{j(\frac{\omega}{2} - \pi)})) $$

Note: This is where the linked paper an therefore image four in the question is confusing. It does not include the linear phase part $e^{-j\frac{\omega}{2}}$ in equation (7). Without this linear phase part, as mentioned in the question, the signal is not a valid factor two down sampled signal. The phase term 'fixes' the equation to make the phase 'repeat' properly according to the new sample frequency.

That's the signal spectrum, and the filter spectrum similarly is:

$$ h[n] \longrightarrow \boxed{z^{1}} \longrightarrow v[n]=h[n+1] \longrightarrow \boxed{\downarrow2} \longrightarrow w[n]=h[2n+1] $$

$$ H_1(e^{j\omega}) = \frac{1}{2} \sum_{k=0}^{1} H( e^{j( \frac{\omega-2\pi k}{2} )  }) e^{j( \frac{\omega-2\pi k}{2} )} $$

In this case $$ H_1(e^{j\omega}) = \frac{1}{2} e^{j\frac{\omega}{2}} \cdot(H(e^{j\frac{\omega}{2}}) - H(e^{j(\frac{\omega}{2} - \pi)})) $$

Note: The linear phase term has opposite sign to the phase term in $X_1$ so these terms will cancel each other when we multiply $X_1$ and $H_1$.

Hence the spectrum at $B_2$ is:

$$B_2(e^{j\omega}) = H_1(e^{j\omega}) X_1(e^{j\omega}) $$
$$ B_2(e^{j\omega}) = \left( \frac{1}{2} \sum_{k=0}^{1} H( e^{j( \frac{\omega-2\pi k}{2} )  }) e^{j( \frac{\omega-2\pi k}{2} )} \right)  \left( \frac{1}{2} \sum_{k=0}^{1} X( e^{j( \frac{\omega-2\pi k}{2} )  }) e^{-j( \frac{\omega-2\pi k}{2} )} \right) $$

in this case:

$$ B_2(e^{j\omega}) = \frac{1}{4} \cdot (X(e^{j\frac{\omega}{2}})H(e^{j\frac{\omega}{2}}) - X(e^{j\frac{\omega}{2}})H(e^{j(\frac{\omega}{2} - \pi)}) -X(e^{j(\frac{\omega}{2} - \pi)})H(e^{j\frac{\omega}{2}}) + X(e^{j(\frac{\omega}{2} - \pi)})H(e^{j(\frac{\omega}{2} - \pi)})) $$

For $A_2$:

$$ A_2(e^{j\omega}) = \frac{1}{4} \cdot (X(e^{j\frac{\omega}{2}})H(e^{j\frac{\omega}{2}}) + X(e^{j\frac{\omega}{2}})H(e^{j(\frac{\omega}{2} - \pi)}) +X(e^{j(\frac{\omega}{2} - \pi)})H(e^{j\frac{\omega}{2}}) + X(e^{j(\frac{\omega}{2} - \pi)})H(e^{j(\frac{\omega}{2} - \pi)})) $$

Note that the aliasing cancellation happens at the reconstruction part of the system: As can now be seen $A_2 + B_2$ results in one channel, whereas $A_2 - B_2$ will produce the other channel.


Best way to detect steps and slopes in a signal

This problem is heavily related to the notion of piece-wise constant or piece-wise linear approximations. Either to fit the experimental curve at best, or to find the optimal parameters may requires some more details, like: do you know the number of parts? Any idea about the noise? Some literature to start with:

Best piecewise constant approximation of a function of single variable
Models and Algorithms for Optimal Piecewise-Linear Function Approximation
Fitting Piecewise Linear Continuous Functions


How to calculate the phase shift AND time delay of non-periodic signals

A pure time delay could be determined by looking for a peak in the cross correlation. But in your case $f2$ might also have an overall phase offset.

You could try to compute two cross correlations:

$$
\begin{align}
x &= cross(f1,f2) \\
y &= cross(f1,hilbert(f2)) \\
\end{align}
$$

where $hilbert(f2)$ refers to an overall 90° phase shifted version of $f2$. If you combine those two like this

$$ z = \sqrt{x^2 + y^2} $$

you should get something that is independent of the phase shift and shows you a peak at the correct time delay $\Delta t$. The "phase" at that peak, $atan2(y,x)$ should give you the phase offset $\Delta\phi$.

I don't know if such a problem is usually solved this way and I have not tried it myself. But it might work.


Why would you use DCT as oppose to DFT for detecting oscillation in a signal?

Just keepin' it $\Re$eal.

The DFT records the timing offset of a particular sinusoid frequency as the phase of the complex coefficient.  Since the DFT of a real signal is conjugate-symmetric, the corresponding negative frequency coefficient gets the same magnitude but negative phase.  FFT routines frequently skip calculating the negative spectrum for real signals.
This leaves NFFT/2+1 non-redundant complex coefficients, however the DC and Nyquist bins are completely real.  It would be possible to combine the DC and Nyquist bin together, making zero growth in data size (for even length signals anyway).

But this kind of flipping between Complex and Real structure may not be convenient for analysis or computation.

The DCT captures essentially the same information about the real signal, except the timing of the sinusoid is encoded as the sum of cosine waves.
The DCT has the advantage of representing a unitary transform without mixing real and complex numbers.  The DFT is also unitary, but it maps $\mathbb{R^N} \rightarrow \mathbb{C^N}$.


Principle of Mapper in frequency diversity in OFDM

In transmitter the N symbols from the series to parallel convertor are repeated M times at the mapper giving you M sets of N symbols. This is a guess, but I  assume the mapper is spreading the contents of each of these M sets out across the whole transmit frequency band so the repeated symbols are not concentrated in any one region. This way if any there is a null in any particular region of the band only a few of the repeated symbols would be affected. This would give you frequency diversity.

AT the receiver the contents of each of the M sets are weighted and combined in a kind of majority rule detection scheme to determine the estimated N symbols.


resampling from fractional to integer sampling rate

You can use the resample command mentioned as follows to accomplish this:

y=resample(x,1250,2541); % x is the input signal


FMCW Radar using FSK Superimposed LFM ramp

First thing first: Once you start talking about stepped FM waveforms with the goal of achieving LFM-like range resolutions you are introducing a more complex receiver along with other associated disadvantages. Stepped-LFM waveforms require that you process each pulse individually, and thus the receiver must reconfigure itself for every pulse which are at different frequencies. Having said that, you must transmit multiple pulses (think multiple PRIs) in order to reap the range resolution benefits. Contrast this to an LFM system that has the large enough instantaneous bandwidth where it can achieve high range resolution with only a single pulse.
This documentation implies that you use the stepped-FM portions of the waveform to achieve high range resolution and the FSK processing to perform Doppler filtering.
Specifically, with this waveform you can process the strictly stepped-FM portions (the solid black lines) as mentioned before to achieve the range resolution required. You can then pass the FSK portions (dotted lines) through a Doppler filter bank, after appropriate mixing and integration, to yield the target velocity.
As you can see, the receiver for such a system sounds pretty complex! In the past, hardware was not up to par to provide the instantaneous bandwidth and computational power that we have today, so great effort was made to design these complex receivers to achieve the desired performance.
Thankfully, hardware has come a long way even since this datasheet was published and we now have access to high instantaneous bandwidth systems which allow for simpler receiver designs. Modern DAC/ADC and embedded processor technology has enabled multi-GHz bandwith waveform generation  and processing. We can easily generate a multi-hundred-MHz LFM pulse and process it digitally via a matched filter using a rather simple receiver design.


Why in formula of power of periodic signal is $2N+1$?

Isn't $-N$ to $N$ $2N+1$ index in total (you probably missed the 0 here)
$$ P = \lim_{N\to\infty} {\frac{1}{2N+1} \sum_{n=-N}^N} \lvert x[n] \rvert^2 $$


efficient spectrum/fft analysis of sparse signals (found in radar applications)

One solution for this kind of problem is called an Instantaneous Frequency Measurement (IFM) Receiver. The pulses are located with an envelope detector and a threshold. The average phase change per sample is computed across each pulse; this is proportional to frequency of the pulse.


You need to be sure about:

Noise power is bounded
The power of the signal step is higher than 4 times that noise power.

If the noise power is bounded and the change in the signal level is higher than 4 times that noise power (from the pictures it seems to be so, I mean, the step is more than twice the max value of noise), why not take the derivative of the signal? something like:

step = abs(x[n] - x[n-1]);
if (step > threshold)
{
  // Event has happened:
  // 1. Take timer value
  event_time = get_time();
  // 2. Update the state machine
  state = (state == START_STATE) ? STOP_STATE : START_STATE;
}

Well, the get_timer() is specific to your MCU and there are better ways to implement a Finite State Machine, but that is the idea. It is fast and fulfills the requirements, I would not go further.


Which Transformation, Or Similarity Metric, Is Rotation, Shift and Scale Invariant?

Use a methods called Error Level Analysis (ELA). It can solve your problem as many forensic applications use it to detected doctored reigns in digital images.


I doubt that Hilmar will agree (maybe he will) but I will offer this answer as a counter to the one that has already been (prematurely) accepted by the questioner.

Hilmar's answer is wrong.  It is factually wrong and analytically wrong.

    1) Does the DTFT take only infinite input sequences?

  No. There is nothing that prevents you from applying the formula to a finite sequence.

That's sorta true but it must be said that the "formula" is explicitly of an infinite sequence.  However, there is nothing in that formula that says any particular sample must take on any particular value (like non-zero) other than a finite value.  There is nothing that prevents all but a finite subset of samples having zero value.

The DTFT is:

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

That's what the formula is and all of the $x[n]$ have a defined value, even if it's zero.

  However, a finite time sequence has an infinite spectrum, so you can't sample it without some amount of aliasing

No one is sampling it (again).  It has already been sampled.  The question is: What happens if the length of the sequence (which is, in the DTFT, infinite in length) is shortened to a finite length?

And the answer is: It depends on the manner in which the sequence is shortened to a finite length.  Or, if the sequence was never infinite in length and later applied to the DTFT, how are the infinite number of samples outside of the original finite length defined.

But if you do define all of those samples outside of your original finite length to be zero, that is saying:

$$ x[n] = 0 \qquad \text{for } x<0 \text{ and } x\ge N < \infty  $$

Then the DFT is equivalent to:

$$ X(e^{j\omega}) = \sum\limits_{n=0}^{N-1} x[n] \, e^{-j\omega n} $$

And then it starts to look like the DFT.  If $X(e^{j\omega})$ is sampled at $N$ equally-spaced values on the unit circle including $X(e^{j0})$, then it is equivalent to the DFT

$$ X[k] = X(e^{j\omega}) \Big|_{\omega = \frac{2 \pi k}{N}} = \sum\limits_{n=0}^{N-1} x[n] \, e^{-j2 \pi k n/N} $$

But something misleading in Hilmar's answer must be corrected:

Applying a finite-length sequence to the DTFT has nothing to do with aliasing.  Nothing at all.


What are poles and zeroes (with respect to the inputs and outputs of a system)?

A transfer function is defined as the ratio of the transform of output to the transform of input where all initial conditions are zero.