Why do we need to make simulation not only on baseband but also on passband?

So, we'll have to look at two things: Passband Signals and Passband Systems.

Let's start with the general notion that the consideration in domain A is equivalent to  consideration in domain B if there's an invertible mapping between A and B that covers all possible inputs and outputs of the consideration.

Signals

That mapping for the complex baseband and real passband domain for signals is well known; mathematically, it's the Hilbert transform, followed by multiplication with a complex exponential of negative center frequency of the passband.

In the other direction, it's the multiplication with a complex exponentionential of a frequency high enough so the complete baseband is at positive frequencies, and then removal of the imaginary part.

That immediately posts three requirements:

The Hilbert transform must be well-defined
There must be a center frequency for the passband signal
There must be a frequency higher than half the bandwidth of the baseband signal.

2 and 3 are equivalent, as far as I can see – you can only use baseband if your signal has been a passband signal before¹.

About the 1. requirement:

Tough one, I've been trying to generally state a domain of definition for the Hilbert transform, but that's nothing for a lazy Sunday noon. Concretely, and function $f(t)$ that fulfills

$$\left\lvert\int\limits_{-\infty}^{\infty} \frac{f(t)}{t-\tau}\right\rvert < \infty$$

is Hilbert-transformable, and you need to account for the fact that you can use all your Lebesgue intergral magic here (you'll really have to, that thing usually has a singularity, and that is an example often used to motivate Cauchy's principal value theorem for electrical engineers).

Now, if $f\in L^p,\, p>1$, then that's a given; but periodic (and hence, not-energy-)signals do work too; tempered distributions do, and a few other things.

Considering the fact that Parseval's theorem already tells us that, inherently, a signal that's limited in energy in frequency domain (and every band-limited signal has to be) needs to be in $L^2$, we can infer that there is none band-limited passband signal that you cannot equivalently observe in baseband.

Systems

This is the tricky part. We're looking for passband systems that we can't map back and forth from passband to baseband.

¹ – that's true for more than signals with closed interval support; compactness suffices, as far as I can tell, but your mixing function then is no longer a single complex exponential. For example, you could have an infinite real tone comb with frequencies $N\pi,\,N\in\mathbb N$, and your mapping to baseband could be Hilbert followed by multiplication with a complex tone comb $\sum_{m\in\mathbb N} e^{jmt}$, and would still end up with a proper baseband, and a way of converting that back to "strangely defined passband". But that's more about generalizing the notion of pass- and baseband than it goes to answer the question...


The graph of $F_a(t)$ is a rectangle with base length $a$ and height $a^{-1}$. The area of the rectangle is therefore equal to $a \cdot a^{-1}=1$. The Dirac delta distribution is often defined as a certain (see below) limit of a sequence of suitable functions of integral 1 in the limit of vanishing support. This definition is already very intuitive, so that your intuition criterium should be met.

A more rigorous look uses a different definition of the delta distribution through the equation

$$\int_{-\epsilon}^{+\epsilon} f(x) \delta(x) dx = f(0)$$

for any open neighbourhood of radius $\epsilon$ around $0$ on which $f$ is smooth, as well as the condition $\delta(x\neq 0)=0$.

If you plug in your $F_a(x)$ in place of $\delta(x)$ we get

$$\int_{-\epsilon}^{+\epsilon} f(x) F_a(x) dx$$

Now the limit of $F_a$ for $a\to0$ is not to be understood of a limit in the usual function sequence sense. That limit does not exist in terms of a function. Rather, the limit is meant to be taken outside an integral under which $F_a$ appears. So taking the limit happens like

$$\lim_{a\to0}\int_{-\epsilon}^{+\epsilon} f(x) F_a(x) dx$$

In order to understand the limit, we can start at any finite but small enough $a$. If we choose $a<2\epsilon$ and use that $F_a$ and with it the function under the integral vanishes outside of $[-a/2,a/2]$ and offer no contribution to the integral, the boundaries can be adjusted accordingly.

$$\lim_{a\to0}\int_{-a/2}^{+a/2} f(x) F_a(x) dx$$

Within these boundaries $F_a$ is constant and equals $1/a$. So the integral simplifies to

$$\lim_{a\to0} \frac{1}{a}\int_{-a/2}^{+a/2} f(x) dx$$

The requirements for $f$ are that it be smooth on the open interval of our integral. So we can expand $f(x)=f(0)+x \cdot f'(0) + \mathcal{O}(x^2)$ which together with the linearity of the integral results in

$$\lim_{a\to0} \frac{1}{a}\left( f(0) \int_{-a/2}^{+a/2} dx + f'(0) \int_{-a/2}^{+a/2} x dx + \int_{-a/2}^{+a/2} \mathcal{O}(x^2) dx\right)$$

The first integral contributes $a$, the second integral vanishes because of symmetry of the boundaries and the odd function under the integral, and the third integral contributes $\mathcal{O}(a^3)$.

This gives

$$\lim_{a\to0} \frac{1}{a}\left( f(0) a + \mathcal{O}(a^3)\right) = \lim_{a\to0} \left( f(0)  + \mathcal{O}(a^2)\right) = f(0) $$

This result tells us that your $F_a$ works as a $\delta$ distribution if the limit is taken outside the integral. Or in other words, the "limit" of the function sequence is a delta distribution.

Edit: I just noticed I did not go into the second condition of $\delta(x \neq 0) = 0$. To show this for $F_a(x)$, just observe that for any $\hat{x}\neq0$ you can find an $\hat{a}$ so that $F_a(\hat{x})=0$ for all $a<\hat{a}$.

Edit 2: Why is the condition $\delta(x \neq 0) = 0$ necessary at all? This one is a bit subtle. The integral condition should uniquely define the distribution since it must hold for all test functions $f$ from the space of locally smooth functions. The problem is that the test functions are not enough to rule out distributions that are identical to $\delta$ only almost everywhere. For example, if we had a distribution $\bar{\delta}$ that coincides with $\delta$ with exception of $\bar{\delta}(1)=1$, then the integral condition would still be met because the set that supports the deviation is of measure zero. The requirement of vanishing $\delta$ for non-zero arguments excludes these cases. It is however often omitted because it does not present a problem unless you let the delta distribution act on other distributions.


The application of a filter to a signal in the time domain is defined as convolution.  Convolution in the time domain is mathematically equivalent to multiplication in the frequency domain. So if we filter a signal, and then filter the output of the first filter, this is equivalent to multiplying the signal by both filters in the frequency domain. Therefor, this is equivalent to filtering with a third single filter whose frequency response is equal to the product of the first two filters’ frequency responses. Because frequency responses are complex numbers, the resulting magnitude is equal to the product of the two input magnitudes and the resulting phase is the sum of the input phases.  This can be easily extrapolated for more than two filters as well.


An intuitive explanation as to why a matched filter is time reversed?

Picture the transmitted signal as a «signature». You want to find some process that maximize the probability of detecting that signature even when there is noise.
What do you do to find some signature buried in noise? You make a template that exactly match the known signature, and you slide it back or forth in time, noting how much the actual signal deviates from the ideal template in any one spot. The time-shift that gives a «large enough» correspondence is your assumed signal location.
The time-reversal of the template is just about reversing the reversal of that argument in the convolution operation. Think of it as non-reversed correlation if you prefer (leaving out complex numbers right now).
-k


convolve to differentiate black and white colors

You have to make all the values zero in the matrix except the middle row and thus, what you can start thinking is how to get such filter matrix which can give you this result. One of the many filter matrices you can think of, which satisfy your requirements is [0 1 0;0 0 0;0 -1 0].


Pitch detection algorithm for minsample size

Consider trying an upsampled or interpolated ASDF, AMDF, autocorrelation or other similar periodicity estimation algorithm.

There in an information theoretic time versus frequency resolution versus noise trade-off.  At a sample rate of 44100, estimating 440 Hz +-2 Hz might require somewhere in the range of 2 to 6 times 44100/440 samples (to determine the existance of periodicity or pitch) up to 2.5 to 4 times 44100/2 samples (to separate 2 spectral lobes), depending on the S/N ratio, and whether that noise environment requires windowing, or allows sub-sample or sub-bin interpolation in your periodicity estimator and/or FFT.

BTW, "super resolution" for tuners is usually considered to be +-1 cent or better, even though this is below normal human audibility except for beating.


Since you care about the time it takes, you'll want to use an optimized FFT implementation. That would be FFTw or FFTS, realistically. Historically, FFTw is way dominant in software (I mean, even Matlab uses that), but FFTS works really well on ARM, so might be the better choice for your Pi.
Then you'd just write a minimal C++ program that

creates an fftw_plan  / ffts_plan_t
generates random data (std::uniform_real_distribution is your friend in generating real and imaginary parts of numbers)
uses using sc = std::chrono; auto start = sc::steady_clock::now(); as a "start" timer,
fftw_execute/ffts_execute the plan on your random data
auto finished = sc::steady_clock::now();, and auto ms = sc::duration_cast<sc::duration<double>>(finished-start).count()*1000; to get the execution time in milliseconds.

You might want to run the actual plan multiple times (maybe even modifying the data in between) and divide the total time by the number of executions, to reduce the measurement error.

A note on the raspberry Pi: If you're using a modern RPi, don't use 32 bit operating systems, but go for an Aarch64 OS. In my experience, not having the SIMD instructions available in 64 bit mode really really hurts numeric throughput performance. (For an FFT, the better loop performance of aarch64 compilers might also help significantly; haven't benchmarked that.)


How do I stabilize my oscillating system?

First, combine your two variable set of first order differential equations into a single variable second order one.

$$ \frac{d^2y}{dt^2} = c \frac{dx}{dt} = acy + bc $$

$$ \frac{d^2y}{dt^2} - acy = bc $$

Solve:

$$ y(t) = C_1 e^{\sqrt{ac} \cdot t} + C_2 e^{-\sqrt{ac} \cdot t} - \frac{b}{a} $$

$$ x(t) = \int{ ( ay(t) + b ) dt }  $$

$$ x(t) = C_1 \frac{a}{\sqrt{ac}} e^{\sqrt{ac} \cdot t} - C_2 \frac{a}{\sqrt{ac}} e^{-\sqrt{ac} \cdot t} + C_3 $$

Now analyze:

Initial conditions have to be set so that $ C_1 = 0 $ or your system will grow exponentially (assuming $ ac > 0 $).  Also, $ C_3 $ will need to be zero.  If $ ac < 0 $ then the system will oscillate forever since there is no damping term in the second order equation.

$$ \frac{dy}{dt} = \sqrt{ac} \left( C_1 e^{\sqrt{ac} \cdot t} - C_2 e^{-\sqrt{ac} \cdot t} \right) $$

$$ \frac{dx}{dt} = a \left( C_1 e^{\sqrt{ac} \cdot t} + C_2 e^{-\sqrt{ac} \cdot t} \right) $$

Initial Conditions:

$$ y(0) = C_1 + C_2 - \frac{b}{a} $$

$$ \frac{dy}{dt}(0) = \sqrt{ac} \left( C_1 - C_2 \right) $$

$$ x(0) = \frac{a}{\sqrt{ac}} ( C_1  - C_2 ) + C_3 $$

$$ \frac{dx}{dt}(0) = a \left( C_1 + C_2 \right) $$

With $ C_1, C_3 = 0 $:

$$ y(0) = C_2 - \frac{b}{a} $$

$$ \frac{dy}{dt}(0) = - \sqrt{ac} \cdot  C_2  $$

$$ x(0) = -\frac{a}{\sqrt{ac}} \cdot C_2 $$

$$ \frac{dx}{dt}(0) = a  \cdot  C_2 $$

It's been a while, I hope I did it right.

Ced

==================================

Followup:

Yes, you can include an $x$ term in the first equation, or a $y$ term in the second equation, or both.  Then repeat the process I went through above to put your system in this form:

$$ \ddot y + A \dot y + B y = 0 $$

Whatever constant you get on the right side will be part of the particular solution.  Add it in later as I did above.

This is the basics of solving a homogenous second degree linear differential equation:

Assume the solution is of the form $ C e^{rt} $

$$ r^2 C e^{rt} + A r C e^{rt} + B C e^{rt} = 0 $$

$$ \left( r^2  + A r + B \right) \left( C e^{rt} \right) = 0 $$

One, or both, of the factors has to be zero.  $ C e^{rt} $ won't be, so:

$$ r^2  + A r + B = 0 $$

This is a standard quadratic equation.  With real valued coefficients you will either get two real roots, or a conjugate pair of complex roots.

$$ r = - \frac{ A }{2} \pm \sqrt{ \left(\frac{ A }{2} \right)^2 - B } $$

This will give you $r_1$ and $r_2$.  Since it is a linear differential equation, any linear combination of the two solutions will also be a solution.

$$ y(t) = C_1 e^{r_1 t} +  C_2 e^{r_2 t} $$

The real part of the root will determine the growth behavior.  If it is negative the system will converge to an equilibrium point, if it is positive, the system will grow exponentially.  The imaginary part, if there, will determine the frequency of the oscillations.  Therefore, you can calculate the parameters from observations by measuring the rate of decay to get the real part(s), and measuring the frequency to get the imaginary part.  Once you have those, you should be able to solve for your $a,b,c$ and $d,e$ if you include them.

Of course, you could leave your equation set in matrix form and solve it with eigenvalues and eigenvectors, but that is a whole 'nuther discussion.

Ced


Deconvolution of 1D Signals Blurred by a Gaussian Kernel

I've explained it once on StackOverflow.

Your signal can be represented as a vector, and convolution is multiplication with an N-diagonal matrix (where N is the length of the filter). I am assuming for the sake of the answer that the filter is much smaller than the signal

For example:

Your vector/signal is:

V1
V2
...
Vn

Your filter (convolving element) is:

  [b1 b2 b3];

So the matrix is nxn: (Let it be called A):

[b2 b3 0  0  0  0.... 0]
[b1 b2 b3 0  0  0.... 0]
[0  b1 b2 b3 0  0.... 0]
.....
[0  0  0  0  0  0...b2 b3]

Convolution is:

A*v;

And de-convolution is

A^(-1) * ( A) * v;

Obviously, in some cases de-convolution is not possible. These are the cases when you have singular A. Even matrixes that are not singular, but close to being singular, can be problematic, as they will have large numeric error. You can estimate it by computing the condition number of the matrix.

If A has low condition, you can compute the inverse, and apply it on the result.

Now, let's see some examples in Matlab:

First, I've made a function that computes the convolution matrix.

function A = GetConvolutionMatrix(b,numA)
    A = zeros(numA,numA);
    vec = [b  zeros(1,numA-numel(b))];
    for i=1:size(A,1)
        A(i,:) = circshift(vec,[1 i]);
    end
end

Now, let's try to see what happens with different kernels:

    b = [1 1 1];
    A = GetConvolutionMatrix(b,10);
    disp(cond(A));

The condition number is :

 7.8541

This one is problematic, as expected. After averaging, it is hard to get back the original signal.

Now, let's try some milder averaging:

b = [0.1 0.8 0.1];
A = GetConvolutionMatrix(b,10);
disp(cond(A));

The result is:

1.6667

That goes well with our intuition, mild averaging of original signal is much easier to reverse.

We can also see how the inverse matrix looks like:

 figure;imagesc(inv(A));

Here is one line from the matrix:

  0.0003   -0.0026    0.0208   -0.1640    1.2910   -0.1640    0.0208   -0.0026    0.0003   -0.0001

We can see that most of the energy in each line is concentrated in 3-5 coefficients around the center. Therefore, in order to deconvolve, we can simply convolve the signal again with this approximation:

   [0.0208   -0.1640    1.2910   -0.1640    0.0208]

This kernel looks interesting! It is a sharpening operator. Our intuition is correct, sharpening cancels out blur.


How can a system be unstable if $L(j\omega)$ is never exactly $-1$?

You cannot make conclusions about the stability of a system by only considering its transfer function evaluated on the imaginary axis $s=j\omega$. Replacing $s$ by $j\omega$ in the transfer function only makes sense for a stable system, otherwise you get a function of $\omega$ that does not describe the system, but another (stable) system.

Let me explain this by an example. Assume a causal system with a transfer function

$$H(s)=\frac{1}{s+a}\tag{1}$$

If $a>0$, the system is stable and by substituting $s=j\omega$, we get its frequency response

$$H(j\omega)=\frac{1}{j\omega+a}\tag{2}$$

The corresponding impulse response is

$$h(t)=e^{-at}u(t)\tag{3}$$

where $u(t)$ is the unit step response. We have

$$H(s)=\int_{-\infty}^{\infty}h(t)e^{-st}dt,\qquad \text{Re}\{s\}>-a\tag{4}$$

and

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

Now imagine that we decrease the value of $a$. The system remains stable as long as $a>0$, but as soon as $a$ becomes $0$, $H(j\omega)\rightarrow\infty$ for $\omega=0$, and the system becomes unstable. If we further decrease $a$, the system remains unstable, but the expression $(2)$ does not go to infinity anywhere. However, the point is that the expression $(2)$ is not the frequency response of the system anymore if $a<0$. The reason is that the integral $(5)$ does not exist anymore. So it is pointless to try to deduce stability or instability from the frequency response of a system, because an unstable system has no frequency response.

What you get if you replace $s$ by $j\omega$ in the transfer function of an unstable system is the frequency response of another system with the same expression for the transfer function, but with a different region of convergence (ROC). The ROC is the region in the complex $s$-plane for which the integral $(4)$ converges. For the causal system in our example, the ROC was given by $\text{Re}\{s\}>-a$ (the ROC is always a right half-plane for causal systems). For an unstable system with $a<0$, the ROC does not include the $j\omega$-axis, and, consequently, $H(j\omega)$ does not exist. In that case, the expression $(2)$ is the frequency response of a different system with transfer function $H(s)$ given by $(1)$ but with a different ROC given by $\text{Re}\{s\}<-a$. That system is stable but anti-causal, and its impulse response is given by

$$h(t)=-e^{-at}u(-t)\tag{6}$$

which is clearly different from $(3)$.

Concerning your second question, note that $e^{\pm j\pi}=-1$, so a phase shift of $\pi$ for certain frequencies implies a sign change, which turns the negative feedback into a positive feedback. This can be compared to the positive feedback in an audio system where a microphone picks up the loudspeaker signal, which is fed back into the amplifier, sent to the loudspeaker, etc.


Implementing FIR/IIR "alike" filter in R

The function filter() in the signal package does exactly what you are asking for. Note that, by default, it uses ARMA filters (a combination of AR and MA). Thanks to that, you can implement both FIR and IIR filters, because the difference equation for an ARMA model is:

$$y(n) = \sum_{m=0}^{M} a(m)x(n-m)+\sum_{k=1}^{K} b(k)y(n-k)$$

To achieve your FIR filter, note that it is enough to set all $b$ coefficients to $0$. For the IIR filter, you should set $a(0)=1$ and every other value of $a$ to $0$.

Another maybe obvious yet important observation is that you can add the noise $v(n)$ at the end of the script, because the noise is white (at least that's what I understand from your question). Thus, there is no correlation between past and present samples of the noise, letting you adding it after applying the filter.


The value of :

img1 = imageio.imread(".jpg")

are rgb value ? i mean are they ranged between 0-255 or between 0 and 1.
Cause if they where all between 0 and 1 then:

((img1 - img2)**2)

will have always value below 10 indeed the maximum value of this operations could be only 1 and therefore the maximum value obtained here :

np.sqrt(squared_diff[:,:,0] + squared_diff[:,:,1] + squared_diff[:,:,2])

could be np.sqrt(3) and that's always minor than 10 and ">10" operator will return an array with all False values


Probability error for antipodal systems and orthogonal systems

For equally likely
binary signals $s_0(t)$ and $s_1(t)$ of energies $E_0$ and $E_1$ respectively
in an AWGN channel with two-sided power spectral density $\frac{N_0}{2}$, the
error probability is given by
$$P_e = Q\left(\sqrt{\frac{E_0 + E_1 - 2 \langle s_0, s_1 \rangle}{2N_0}}\right)
\tag{1}$$
where
$\displaystyle \langle s_0, s_1 \rangle = \int_{-\infty}^\infty s_0(t)s_1(t)\,\mathrm dt$ is the inner product of the two signals.

For antipodal signals, $E_0 = E_1 = E$ and $\langle s_0, s_1 \rangle = -E$ giving
everyone's favorite answer $Q\left(\sqrt{2E/N_0}\right)$ while for equal
energy orthogonal signals, the inner product is $0$ giving everybody's second
favorite answer $Q\left(\sqrt{E/N_0}\right)$. For some details of the proof
of $(1)$,
see this answer that I wrote
some time ago.  For even more detail,
see
this Lecture Note of mine.


1) A Fourier serie is written as $$\frac {a_0} 2 + \sum \limits _{m=1} ^\infty (a_m \cos \frac {2 \pi m t} T + \mathbf{i} \, b_m \sin \frac {2 \pi m t} T)$$
Unless your signal is even, i.e. $f(-t)=f(t)$ the fourier decomposition will have a real (cosine) part and an imaginary (sine) part. The middle bottom plot is more likely to be the power spectral density i.e. the fourrier transform's squared modulus.

2) Indeed, the Fourrier basis is a orthonormal basis, in a mathematical sense.
$$
\zeta _m(t)=\cos \frac {2 \pi m t} T + i \, \sin \frac {2 \pi m t} T
$$
$$
\int_{0}^{T}{\zeta _k(t)\zeta _l(t) \, dt}=\delta _{kl} %\quad \forall k,l\in[0,T-1]
$$

The two conclusions above still holds in the 2D case simply use as basis functions:
$$
\zeta _{u,v}(x,y)=e^{-i2\pi(ux+vy)}
$$
Define an image $f(x,y)$ and it's Fourier transform :
$$
F(u,v)=\int \int f(x,y) \zeta _{u,v} \, d x d y
$$
$F(u,v)$ isn't real unless $f(x,y)$ is even respect to $x$ and even respect to $y$, and so when building the Fourier transform's representatio, it is more likely that the represented value is the Fourier transform's squared modulus.
In the 2D case, Fourier's basis function still form a basis, it is complicated to define a base in a functional space. The two main things to note are :
First, the scalar product between basis functions $\zeta _{u_1,v_1}$ and $\zeta _{u_2,v_2}$ is always equal to 0 unless $u_1=u_2$ and $v1=v_2$ then it is equal to 1 (orthonormality). One must then prove that every function $f(x,y)$ is uniquely writtable as a linear combinaison of $\zeta _{u,v}$'s. It is true for periodic function, but I can't provide you the demonstration.

I let you check the definition of a basis for a vector space, if you understand it it will give you a good grasp of what is a basis in general.


How to filter a signal using a bandpass filter consisted of two moving average filters?

The Japanese link actually implies how to derive Eq.5, and ironically they also refer to an existing dsp.se answer at the bottom.
Derivation of Eq.5 is as follows:
Consider a moving average filter of length $N$, with the impulse response $h[n]$:
$$h[n] = \begin{cases} ~~~1/N~~~,~~~n=0,1,...,N-1 \\ ~~~0~~~,~~~ \text{otherwise} \end{cases} \tag{1}$$
Magnitude of its frequency response $~H(\omega)~$  (DTFT of $h[n]$) can be shown to be:
$$ |H(\omega)| = \frac{1}{N} \left|\frac{ \sin(\frac{\omega}{2}N)}{\sin(\frac{\omega}{2})}\right| \tag{2}$$
Now, after replacing the $\sin()$ functions with their Taylor expansions in powers of $\omega$, it can be shown by polynomial long division that Eq.2 is also given by :
$$|H(\omega)| = 1 + \frac{1}{24}(1-N^2) \omega^2 + H.O.T. \tag{3}$$ where H.O.T. refers to higher order terms in powers of $\omega$. An approximation for small values of $\omega$ is obtained by neglecting H.O.T.:
$$|H(\omega)| \approx 1 + \frac{1}{24}(1-N^2) \omega^2 \tag{4}$$
Using this approximate frequency response magnitude, we can obtain an approximate cutoff frequency $\omega_c$ at which the magnitude $|H(\omega_c)|$ falls to $1/\sqrt{2}$ of its value at $\omega = 0$, which is $H(0) = 1$.
$$ |H(\omega_c)| = \frac{1}{\sqrt{2}} = 1 + \frac{1}{24}(1-N^2) \omega_c^2 \tag{5}$$
Replace the discrete-time frequency $w_c$ by $w_c = 2\pi f_c /f_s$, where $f_c$ is the analog cutoff frequency in Hz, and $f_s$ is the sampling frequency in Hz. Finally solving the resulting algebraic expression for $f_c$ yields the formula that you refer to as Eq.5 in the document :
$$ |H(\omega_c)| = \frac{1}{\sqrt{2}} = 1 + \frac{1}{24}(1-N^2) \left( 2\pi \frac{f_c}{f_s} \right)^2 \tag{6}\\\\$$
$$ f_c = \frac{1}{\pi} \frac{\sqrt{6 - 3\sqrt{2}}}{\sqrt{N^2-1}} ~f_s ~~ =~~  \frac{0.422}{\sqrt{N^2-1}} ~f_s \tag{7} $$
Eq.7 above is the formula that provides an approximate cutoff frequency calculation for the moving average filter of length $N$ (order $N-1$). In your posted link, Eq.5 there's a slight variation at the scale of $0.442$ instead of $0.422$, probably they tried to apply some correction to the actual cutoff vs approximated one.
Note that in the derivation, we've used an approximation of the DTFT magnitude which was valid as long as $\omega$ was small compared to $\pi$. This means that the approximation will be satisfactory if $\omega_c$ is close to $0$, or in other words, $f_c$ is a small compared to $f_s$. And indeed this will be the case for high order moving average filters. And the approximation gets better as $N$ increases.
Using such two moving average filters with approximate cutoff frequencies $f_{c1}$ and $f_{c2}$ to create a bandpass filter will not be very satisfactory unless your out of band signal energy is insignificand after 20 to 30 dB of attenuation.


How to get sensor covariance matrix from CSI values?

The Channel State Information (CSI) matrix is calculated over a packet of known signals that is sent as part of the communications protocol.

The CSI matrix contains two things (per carrier). The bearing of the "target" and the "strength" of transmission along the line of that bearing. So basically, this can tell you where is the transmitter with respect to the receiver array.

The CSI matrix is calculated via cross correlation between the signals received by each antenna and correlation is normalised covariance.

So, strictly speaking, it is impossible to derive covariance from the CSI because you need access to the received signals themselves. But, the angle (bearing) component of the CSI could be considered as being proportional to the real covariance for the pilot tones that are used to generate them. You could even obtain some theoretical covariance from the pilot tones by delaying them progressively but that would not include the covariance introduced due to noise.

For more information please see here and here.

Hope this helps


Gaussian Blur In MATLAB and Connection to Image Resolution

When applying Gaussian Blur on an image you should care for the STD only.
The rest should be set as following:

STD_TO_RADIUS_FACTOR = 5;
kernelRadius = ceil(STD_TO_RADIUS_FACTOR * kernelStd);
kernelLength = (2 * kernelRadius ) + 1;
mGaussianKernel = fspecial('gaussian', [kernelLength, kernelLength], kernelRadius);

Now, regarding the resolution of the image, this is tricky since it has many parameters in it.
If you want to measure it in Cycles in Spatial (As DFT), you can see the spectrum properties (Cut Off Frequency) of a Gaussian Kernel.

Update

The STD_TO_RADIUS_FACTOR should be determined by the balance between accuracy and performance.
The above is basically a truncated version of the Gaussian Kernel (Which is infinite). The balance sets how good is the approximation vs. how long the Convolution Kernel is (The longer the higher computational resources it requires).


Is there any modern information about a 12-hour periodicity in sonar noise?

Some further googling seems to have shed some light: it appears that shipping traffic is highest at these hours, as demonstrated by Figure 3 in
http://cetus.ucsd.edu/Publications/Publications/McKennaIslandsSymp2009.pdf.
Why that is, I'm not sure. Also, this publication isn't about global shipping patterns, but only that in a certain region on the west coast of the USA - but it certainly goes some way to explaining the noise detected by sonar.
I know this was a very niche question, but hopefully this is helpful to someone who may stumble on my question in the future.


Problem understanding the Expectation Operator

Consider a discrete-time, discrete valued random process $X[n,s)$ whose samples $x[n]$ at integer index $n$ belong to range space $R_{X_n} = \{x_1,x_2,...,x_n\}$ of discrete valued random variables $X_n$ indexed by $n$, such that at each $n$, the value $x[n]$ is the result of the mapping by $X_n$ into $R_{X_n}$; i.e., $X[n,s) = X_n(s) = x[n]$, where $s$ is the sample output of the experiment performed at time $n$.

The expectetation of a particular random variable $X_n$ of the random process is
$$\mathcal{E} \{ X_n \} = \sum_{x_k \in R_{X_n}} x_k P_{X_n}(x_k) = \int_{-\infty}^{\infty} x f_{X_n}(x) dx $$
where those $x_k$ belong to the range space of the random variable $X_n$ and $P_{X_n}(x_k)$ are the probabilities of those values. Note that the integral definition is also available for those who would still like to use a continuous random variable approach with impulsive PDFs; $$f_{X_n}(x) = \sum_k P_{X_n}(x_k) \delta(x-x_k) $$

The correlation between two (real) random variables $X_n$ and $X_m$ is defined as
$$\mathcal{E} \{ X_n X_m \} = \sum_{k} \sum_{l} x_k x_l P_{X_n,X_m}(x_k,x_l) $$
where $P_{X_n,X_m}(x_k,x_l)$ is the joint probability mass function of the two random variables, and signifies the probability of the pair $(x_k,x_l)$ in the range space of the joint mapping.

When those two random variables belong to the same random process, then the correlation is defined as auto-correlation of the random process.

Note that the second sum in your post does not represent an expectation but rather a possibly ergodic estimate of it for a WSS random process.


It depends on the type of noise and type of signal.  Show an example if you want a good answer.  But, that said, in general you probably want to low-pass filter the signal.  If I were you, I'd take a Fourier power spectrum to see if most of the noise is high frequency, and the signal I care about mostly in a lower range.  If they overlap, oh well that's life. I would have to think more about things.

One low-pass filter that's good for noisy signal in many cases is the Savitzky-Golay filter.  It is described in Numerical Recipes, and for Python there's a function in the Python Numpy Cookbook.  It is merely a convolution with a small kernel.   You pick the window size based on the width of the peaks or other features, wide enough to mush out the noise, but not wider than the features.  It can be small, say five points, or bigger like dozens, a hundred maybe.

You also pick a polynomial order - usually I use 2 or 4.  Order 2 is fine for when the window is small, < 10 points or it spans less than half a cycle or so (if your signal resembles a sine) while order 4 is better at matching distorted peak shapes, but likes to have around 9 or more points.  But a lot depends on the shape and frequency of the noise.

As other say in the comments, finding derivatives probably isn't the best strategy, but if you want find derivatives anyway, the Savitzky-Golay filter can do that - simultaneously smoothing and reporting the derivative instead of the signal.


Without more knowledge on $X$, we cannot say a lot. However, if it is IID with mean $\mu$ and variance $\sigma^2$, the sample mean $\hat{\mu} = \frac{1}{N}\left(\sum_{n=0}^{N-1}x_n\right)$ is unbiased, since:
$$E[\hat{\mu}] = E\left(\frac{1}{N}\sum_{n=0}^{N-1}x_n\right) = \frac{1}{N}\sum_{n=0}^{N-1}E\left(x_n\right) = N\mu/N=\mu\,.$$
And its variance goes to zero when $N$ increases:
$$V[\hat{\mu}] = V\left(\frac{1}{N}\sum_{n=0}^{N-1}x_n\right) = \frac{1}{N^2}\sum_{n=0}^{N-1}V\left(x_n\right) = N\sigma^2/N^2=\sigma^2/N\,.$$
Thus, the expectation converges to the actual mean, and the variance of the estimator tends to zero as the number of samples grows. Under these definitions, the sample mean is a consistent estimator. Note that one could try to use other hypotheses: alternative norms, convergence in law, etc.


Edge detection in spatial or in frequency domain?

What you claim is simply not true. One step in edge detection is often a simple high-pass (or band-pass) filter,  and those are commonly applied using fast convolution, which is a frequency-domain algorithm.


How to calculate $p_c$ (binomial probability)

You can use the binomial approximation:
$$
(1 + x)^\alpha \approx 1 + \alpha x,
$$
which should be a very good approximation. The conditions for this being accurate are that $|\alpha x| \ll 1$.

When you substitute this expression into your total probability, you now have a polynomial you can solve for the roots of.


Is there any Spread Spectrum system that uses Rake with cyclic Prefix?

you can read the following paper, they are using almost the same idea, the title of the paper: Symbol Cyclic Shift Equalization PAM-OFDM - A Low complexity CP-free OFDM Scheme ..  Good luck


How to understand the discrete and continuous convolution in matlab?

As Hilmar said you need to be careful of the difference between circular convolution and linear convolution. Linear convolution in the time domain translates to multiplying discrete time fourier transforms (DTFT) or continuous time fourier transforms (CTFT) in frequency and I'm sure you've seen this done.

But as you know, the DTFT and CTFT are infinite sums/integrals so cannot actually be computed by a computer. The discrete fourier transform (DFT) can be computed by MATLAB though and it uses an efficient algorithm called the fast fourier transform (FFT) which runs on the order of NlogN time, where N is the length of the signal. You can look up that multiplying two DFTs is equivalent to performing circular convolution which just "wraps" the signal around (read Oppenheim's DSP book or wikipedia).

Check this out which shows the difference between performing linear convolution, how to use the fft() function in MATLAB to get the result, and how circular convolution is different.

clear,close all;

a = [zeros(40,1); ones(20,1); zeros(41,1)];
b = [zeros(50,1); ones(30,1); zeros(21,1)]/2;
c = conv(a,b);
nfft = length(a)+length(b)-1;
A = fft(a,nfft);
B = fft(b,nfft);
C = A.*B;
d = ifft(C);
e = ifft(fft(a).*fft(b));

figure
subplot(6,1,1)
plot(a)
title('a[n]')
subplot(6,1,2)
plot(b)
title('b[n]')
subplot(6,1,3)
plot(c)
title('Convolution in time domain')
subplot(6,1,4)
plot(d)
title('Convolution by multiplying in freq. domain')
subplot(6,1,5)
plot(e)
title('Circular convolution by multiplying in freq. domain')


Find the frequency response of the system

OK, in the comments you've shown that you can actually do it by yourself. This answer just summarizes the steps and shows you how to write your (correct) result in the form that you were given as an answer.

First of all, note that you have two identical systems in series. If $G(\omega)$ is the frequency response of the system up to and including the first integrator (from the left), then the total frequency response is

$$H(\omega)=G^2(\omega)\tag{1}$$

The response $G(\omega)$ can be seen as a multiplication of two subsystems, the system up to and including the adder, and the integrator. The first of the two has a frequency response

$$G_1(\omega)=1-e^{-j\omega T}\tag{2}$$

Since $G_1(\omega)$ has a zero at $\omega=0$, we don't need to bother with any problems of the integrator at $\omega=0$ (I'm talking about the delta impulse). So for calculating $G(\omega)=G_1(\omega)G_2(\omega)$, we can simply use

$$G_2=\frac{1}{j\omega}\tag{3}$$

which gives

$$G(\omega)=\frac{1-e^{-j\omega T}}{j\omega}\tag{4}$$

and, from (1),

$$H(\omega)=\frac{(1-e^{-j\omega T})^2}{(j\omega)^2}\tag{5}$$

In order to get from (5) the answer that you were given, use this trick to rewrite $G_1(\omega)$:

$$G_1(\omega)=e^{-j\omega T/2}(e^{j\omega T/2}-e^{-j\omega T/2})=2je^{-j\omega T/2}\sin(\omega T/2)\tag{6}$$

Using (6) you can rewrite $G(\omega)$ as

$$G(\omega)=2e^{-j\omega T/2}\frac{\sin(\omega T/2)}{\omega}=
Te^{-j\omega T/2}\frac{\sin(\omega T/2)}{\omega T/2}\tag{7}$$

Squaring (7), using $\omega=2\pi f$ and $\text{sinc}(x)=\sin(\pi x)/(\pi x)$ will give you the final form of the result.


Examples of Causal and Complex Wavelets?

You can always just build your own, for instance by taking the Morlet wavelet shift it to the right until, say, 99% of the area below the function is on the positive axis and then set everything on the negative axis to zero. But bear in mind that it's impossible to build a wavelet (or filter which is the same) that is causal and analytic, i.e., a filter that erases all negative frequencies, which you actually need to get that property you need (separation of phase and amplitude). However a pretty good approximation is always possible.


How to Estimate the Noise of an Image / Estimation of the Noise in an Image?

There are various methods to estimate the noise of a signal (or adapt a threshold to be used later in PCA for example).

There are adaptive methods and statistical methods. Statistical methods employ a statistical model of the noise and estimate from the data, while adaptive methods iteratively filter the data until a certain threshold of (reduced signal) accuracy is reached. These are just generalities of the most used approaches.

The following references deal (mostly) with images from cameras:

Noise Estimation from a Single Image
Automatic Estimation and Removal of Noise from a Single Image
A noise-estimation algorithm for highly non-stationary environments


This may be caused by "border effects" related to subsampling and signal length with perfect reconstruction filter banks.

For 2-band wavelets, there are standard  non-expansive solutions (for instance with periodic border extension) when $2^N$ divides the length $L$ of your signal. In your case, one or two levels work without length increase, becauses $20=20/2+20/2 = (10/2+10/2)+10$. Then, starting at level $3$, you cannot divide $5$ into equal parts, since $2^3 = 8$ does not divide 20.

So:

It signifies that you have decomposed too deep to keep the number of samples constant.
The increase happens because of your signal length, your decomposition  level and your border extention.

For some applications (outside compression) this is not really an issue.


The YIN Algorithm could be of interest to you. It's based on autocorrelation and is used for estimating the fundamental frequency.

A. de Cheveigné, H. Kawahara: YIN, a fundamental frequency estimator for speech and music, Proceedings of ICMC 2005. Download PDF.

A C++ implementation is also available.


Compute time constant of discrete exponential decay pulses

Take four sample values.

$$ y_1 = x[p] $$
$$ y_2 = x[q] $$
$$ y_3 = x[p+d] $$
$$ y_4 = x[q+d] $$

Subtract the third from the first.  Notice that the $C$s cancel out.

$$ y_1 - y_3 = x[p] - x[p+d] = Ae^{-p/\tau}-Ae^{-(p+d)/\tau} =  Ae^{-p/\tau} \left( 1 - e^{-d/\tau} \right) $$

Similarly, subtract the fourth from the second.

$$ y_2 - y_4 = Ae^{-q/\tau} \left( 1 - e^{-d/\tau} \right) $$

Take the quotient of these two differences.  Notice that the $A$s cancel out.

$$ \frac{y_1 - y_3}{y_2 - y_4} = e^{(q-p)/\tau} $$

Take the log of each side.

$$ \ln \left( \frac{y_1 - y_3}{y_2 - y_4} \right) = \frac{q-p}{\tau} $$

Solve for $\tau$.

$$ \tau = \frac{q-p}{ \ln \left( \frac{y_1 - y_3}{y_2 - y_4} \right) } $$

I did this more thoroughly in another answer, but this should give you start.  You want the differences to be large compared to the noise level, meaning d should be large.  $(q-p)$ being larger will help too.  I am guessing that four evenly spaced samples are as good as anything.  Something to investigate.

If noise is a problem you can spread the calculation across more points.  You can either do an average of neighbor points for each $y$ value as long as you use the same weighting or employ more points directly in the quotient as in my referenced answer.

The other answer will also show you how to find the best fit values of $A$ and $C$.

Using more points to mitigate noise.

Here is a way you can fully utilize all your samples to reduce the effects of noise I believe as much as possible.

Select the number of samples ($N$) to be a multiple of four.  Partition the samples into four quarters.  Then set the $y$ values to be the sum of the samples of the corresponding quarters.  You don't need to take an average (an extra division) because that washes out in the taking of the quotient.

$$ q-p = N/4 $$

In addition to taking one measure of the entire curve, you might want to try a sliding window taking many smaller measures.  The values you get should be fairly constant.  If not, your assumption of exponential decay may not be correct.

Sample Python code:

import numpy as np

#================================================
def main():

#---- Set Parameters

        A   = 1.234
        C   = 5.678
        tau = 24.08016032

        N = 100

#---- Construct the pulse

        x = np.zeros( N )

        for n in range( N ):
            x[n] = A * np.exp( -n / tau ) + C

        print x

#---- [Add noise here]

#---- Example Using a Single Point

        p = 10
        q = 20
        d = 20

        y1 = x[p]
        y2 = x[q]
        y3 = x[p+d]
        y4 = x[q+d]

        S = ( y1 - y3 ) / ( y2 - y4 )

        #++[Check S here]

        tau_calc1 = ( q - p ) / np.log( S )

        print tau_calc1

#---- Example by Quarters

        y1 = np.sum( x[p:p+10] )
        y2 = np.sum( x[q:q+10] )
        y3 = np.sum( x[p+d:p+d+10] )
        y4 = np.sum( x[q+d:q+d+10] )

        S = ( y1 - y3 ) / ( y2 - y4 )

        #++[Check S here]

        tau_calc2 = ( q - p ) / np.log( S )

        print tau_calc2

#================================================
main()


Both Cooley-Tukey and Radix DIT & DIF are based on the same principle, dividing the N samples into two groups, and doing the same for the resulting two groups recursively. DIT and DIF generally use Radix2, that is, split N into two N/2 groups and provide a $N\log{N}$ time, while Cooley-Tukey is a generalization which splits it into $N_1$ and $N_2$ groups, where $N=N_1\cdot N_2$. Cooley-Tukey also provides $N\log{N}$ time in most cases.

Now, Radix DIT and DIF (both special cases of the Cooley-Tukey) are very similar and provide the same optimization in terms of operation savings, with the difference being bit reversed samples either on the input or the output. These can be done in $O(N)$ time, either before (DIT) or after (DIF) calculating the DFT.

Moreover, in some cases you don't even have do to the bit-reversal. Some applications (such as convolution) work equally well on bit-reversed data, so one can perform forward transforms, processing, and then inverse transforms all without bit reversal to produce final results in the natural order.

Finally, as an answer to your question, you can use Radix2 DIT or DIF FFT, depending on the application (see above paragraph), with the same result: Instead of:

$4N^2 \quad\text{real multiplications}\\
4N^2-2N \quad\text{real additions}.$

We achieved:

$2N\log_2{N}-7N+12 \quad\text{real multiplications}\\
3N\log_2{N}-3N+4 \quad  \text{real additions}.$

Please note that Radix DIT and DIF algorithms require that $N=2^p$, which is generally easy to satisfy. Increasing the radix gives us $\log_4$ for radix4, etc., which increases the calculation complexity on hardware, and in turn decreases the number of operations by ~$25\%$. Radix 2 and 4 are considered the most common, while Radix 8 (~$8\%$ optimization) and up generally demand too complex hardware for too small optimizations. Cooley-Tukey gives similar results, with more room for tweaking - in case of a specific $N$, for example.

For further optimizations, there are the split-radix algorithm which uses both Radix 2 and Radix 4 , and many more application-specific algorithms.

Update: Other algorithms include:

Bruun's algorithm:

It's not used much because Cooley-Tukey algorithm offers at least the same, and often better optimization. It was developed as a means to calculate DFT of real data (as opposed to complex).

Rader's algorithm:

Rader's algorithm is used for cases where $N$ is a prime number, often a big prime number. Note that above mentioned algorithms require $N = 2^k$, so they don't really work for other $N$. Rader's algorithm guarantees $N\log{N}$ time for prime $N$.

Bluestein's algorithm:

Bluesten's algorithm calculates the DFT via convolution, and is another way of dealing with $N \neq 2^k$. For prime numbers, it offers $N\log{N}$ time, for even numbers Cooley-Tukey algorithm is better, and for odd numbers it offers a somewhat better time than Cooley-Tukey.

Please note, however, that for all such cases for $N$ it's possible to zero-pad the input series until we get an $N = 2^k$ and then use the Cooley-Tukey's algorithm. Of course, for prime $N$ it's more efficient to use Rader's or Bluestein's algorithms.

Moreover, these algorithms were developed at different points in time, so some are more efficient than the others.


Camera Calibration using Single Input Image

Take a look at this paper here: "Straight lines have to be straight" by Faugeras et al.
http://link.springer.com/article/10.1007/PL00013269#page-1

It is straight forward to implement, but essentially removes distortions by observing the straightness of lines.

It should be straight forward to implement if you remove some automation (e.g. you click on the lines manually rather than relying on image processing).


How to implement Matlab welch in python with the same result?

As written in Matlab pwelch document, by using an integer for window, a Hamming window of length window(integer) will be used. Also the nfft in python should set to the next power of two greater than length of window.

freqs, PSD = signal.welch(data, fs=2000, window=signal.windows.hamming(len(window)), nfft=num_fft, scaling='density', return_onesided=True, detrend=False)


Structuring Kalman filter for tracking problem where only position is known

Just a slight nit $\sigma^2$ is a variance $\sigma$ is a standard deviation.

Let me answer your last question first.  In a KF, one has measurements and states, and they are not usually the same.  In a Savasky Golay filter, one is estimating the derivatives of the measurement, not the state.  In your case the measurements are trivially related to the states.

I'm guessing but I think you want something like
$$
\left| \begin{array}{c} x \\ y \\v_x \\ v_y \end{array}\right| =
\left| \begin{array}{cc}
1 & 0 & \Delta_t & 0 \\
0 & 1 & 0 &\Delta_t \\
0 &  0 & 0 & 0 \\
0 & 0 & 0 & 0 \end{array} \right| +
\left| \begin{array}{c} 0 \\ 0 \\ w_1 \\ w_2 \end{array} \right|
$$
for the state model  $\Delta_t$ is your time step and $w_1, w_2$ iid Normal noise.

If this model is your problem,  yes the KF is what is wanted.

What you will need to figure out, is what is process noise and what is measurement noise.

A KF needs measurement noise to work, process noise is optional, but you would need a different state model if you want to just use measurement noise.


please give the reason why every notation for DFT is valid?

It simply does not matter and it is just a matter of convention which version is used. I always use the first definition, which is also the one you will find in (almost) all textbooks on DSP. If you take always the same time domain signal $x[n]$ then it is obvious that $X_2[k]$ (according to notation 2) and $X_3[k]$ (according to notation 3) are simply scaled versions of $X[k]$ (according to notation 1). The one with scaling $1/\sqrt{N}$ is the unitary version of the DFT (i.e. the DFT becomes a unitary transformation). It's just notation number 4 that is a bit weird. In this case, $X_4[k]=X[-k]$, but I have never seen the sign convention like this. In this case, $x[n]$ can be interpreted as the Fourier series coefficients of $X[k]$. But again, all four versions are OK and you're free to use any of them as long as everybody knows which version of the DFT you mean.


Obtain the complex envelope of a two-tone signal

You can represent any bandpass signal by its complex envelope.

You must, however, pick a frequency for which you do that. In your first example, that frequency is $f_c$; in your second example, you need to pick $f_1$, $f_2$ or any other frequency as $f_c$ – and then represent your cosines into products of $e^{j2\pi (f_1-f_c) t}$ and $e^{j2\pi f_c t}$.

Any textbook that introduces the concept of complex envelope (and, usually, complex baseband) should introduce that concept; if yours doesn't, it might be worth having a second textbook for cross-literature.


Yes, that's correct.

The impulse response of a cascade of two systems is the convolution of the individual Impulses responses. Convolution is commutative, so it doesn't matter in which order you convolve.

Impulse responses are ONLY defined for linear time invariant systems, so this will work for systems like equalizers or filters but NOT for things like dynamic processors or limiters.


I would really recommend that you compute the DTFT of the sequences $x_1[n]=a^nu[n]$ and $x_2[n]=-a^nu[-n-1]$. This is very straightforward if you use the formula for the geometric series:
$$\sum_{n=0}^{\infty}q^n=\frac{1}{1-q},\qquad |q|<1\tag{1}$$
It will be important to consider the magnitude of $a$. Note that it's not for nothing that the table in your question specifies a condition on $|a|$.


What is the frequency response function for the Bode Plot pictured here?

Your Bode plot is a straight line. As we are on a semilog-scale, the equation is given by

$$y = m \log(x) + b$$

In this case, $y$ would be the magnitude of the frequency response $|H(j\omega)|$ in $\mathrm{dB}$, and $x$ would be the frequency $\omega$. Taking two points of the plot and solving the $2\times2$ equation system should give you the function you're looking for.

As you didn't put the values of the $\omega$-axis, I'm going to assume that the plot starts at $\omega=1$. In this case, when $\omega = 1$ the magnitude is $0 \ \mathrm{dB}$, so:

$$0 = m \log(1) + b\implies b=0$$

Also, when $\omega = 10$, the magnitude is $20 \ \mathrm{dB}$:

$$20 = m \log(10) +0\implies m=20$$

So the function is:

$$|H(j\omega)|_{\mathrm{dB}}=20\log(\omega)$$

Note that if your plot doesn't begin at $\omega=1$ (what I just assumed to be true), the reasoning would be the same but $b$ would not be zero.


The Zero-Crossing rate threshold for a voiced/unvoiced decision

Plain ZCR criterion is not enough for robust and accurate voiced/unvoiced separation. That being said, your threshold should be adaptive, there is no fixed threshold that works well for all speech waveforms. It depends on the approach you follow but a statistical threshold should do the work most of the time. You can google your question to find relative papers and see if their approach suits your purpose.

Just a note, ZCR is not the best choice for pitch estimation - it just gives a rough estimation.


If sampling rate is high does that increase the error or decrease the error?

Well, i don't see how is that related to DSP or DIP, except of course the issue of sampling, but here is what is going on:

If you are tracking the mouse by latching on an OS' provided "onMouseMove" signal then you are already enjoying (almost) the maximum sampling resolution offered by the combination of your (hardware + OS). That would be the absolute max in Fs accessible by other programs running on the system.

Is that enough? It is, if you are free drawing with small movements trying to get the contour of an object in some photo editing software (for example). But if you click on the mouse and go crazy with it you might start to notice the trick behind the "smooth lines" which is to compose them of linear segments. Therefore, because of the large actual movement compared to the Fs you will start noticing some long lines appearing as the linear segments connect points that are far appart.

Can this Fs be changed? (Increased / Decreased) Yes, it probably can but that would depend on your hardware + software combination. For example, here are a few pieces of information regarding the PS2 mice where you can see a specific command that can be sent to it to change the sampling resolution (search for "set mouse sample rate"). How you actually get to send this command to the mouse so that the sample rate is changed across the system, depends on your OS).

Will higher Fs lead to improved tracking of the actual mouse movement (smaller error)? Yes, that is the point of increasing the Fs, you can track the continuous physical phenomenon much better.

Why don't we crank the Fs up all the way then? Because depending on your Hardware + OS there will be a limit that you can increase the Fs before it starts becoming "noticable" by your system. What you want to avoid is having the OS triggering the "onMouseMove" signal / message so fast that it floods the message queue and at the extremes the system ends up absorbed to serving just this kind of requests.

Similarly, if you are tracking the mouse position with any other means (a camera looking at the screen for "an arrow that points up left" for example, or the same thing by screen grabbing) and possibly with another $Fs_{Track} < Fs_{Sys}$ then the same thing applies. However, in this case, MAYBE you can somehow increase $Fs_{Track}$ by improving the scanning process...but the best case scenario is that you will be approaching $Fs_{Sys}$ not exceeding it.

I hope this helps.

P.S. Regarding "every flicker", modern mice will report a resolution in DPI (Older mice had mechanical rotary encoders and just good enough resolution). Therefore, you can consider your mouse reporting its position over a very fine canvas. In general, mice are good in discriminating movement and will report robust positions over this fine canvas. There are exceptions however where an optical mouse, at rest, will report changes in the position depending on the surface it rests on. In this case, the mouse gets "tricked" and you can avoid this flicker by changing the mat it rests on. For any other means of tracking, the "flicker" will depend on the tracker and Signal to Noise Ratio (SNR).


Recommended resources for noise reduction

Signal denosing is a well-studied technique in signal processing. It first began using simple techniques such as filtering. In this approach, the emphasis is laid on designing filters which can perform denosing techniques in a fast and efficient manner.

Later, when wavelet theory was developed, some researchers used wavelets to denosing the signal. Essentially, we can apply thresholding techniques (soft or hard threshold) on the coefficients of the wavelets. On reconstruction using the dominant coefficients, one can denoise the signal.

In the early 2000s, sparsity-based techniques were introduced where the input signal is represented using an over complete dictionary. In the reconstruction step, sparsity inducing norms were used to extract only the dominant columns of the dictionary. Popular approaches include LASSO and total variation denosing.

Reference

Useful reference with MATLAB codes http://www.numerical-tours.com/matlab/


Given a wide band signal with equal SNR over the whole spectrum, what should the stop band attenuation be for a decimation filter?

Imagine an incoming signal with 60dB SNR over the whole spectrum.

The 60 dB isn't what matters here – what matters is that you've got the same power spectral density all over your spectrum, including the 39/40 that you'll alias onto your remaining band.
So, if you don't want the SNR of that remaining band to be affected, you'd need infinite attenuation; can't have that.
You'd need to define the amount of SNR degradation you'd want to accept. If you didn't filter at all (anti-aliasing filter with a 0 dB stopband attenuation), you'd get -16 dB (that's roughly 1/40) of SNR. The fact that you had -60 dB of noise immediately stops mattering, because you're aliasing 39 times as much power onto your band as your band has.
If you want to, say, don't reduce your SNR to below 57 dB, well, you need your 16 dB stronger aliasing energy to be suppressed to -60 dB (then you have as much noise in your band as aliases, and that halves your SNR from 60 to 57 dB), that's 76 dB of stopband attenuation.
That is a lot. It's actually hard to implement that reliably and performantly. So, you'd often design a good halfband filter with high suppression and a decimation of 2, and concatenate a couple of these before finally doing the rest with one or again cascaded FIR filters.


The Euclidean Distance Transform can produce dilations and erosions with suitable parameters and filtering.

The algorithm that Photoshop uses, and the one that is best suited for stroking, is to calculate the Euclidean Distance Transform in integers and without taking the square root (i.e. calculate distance squared). This can be made extremely fast using the technique from Meijster.

To get anti-aliasing, just use the fractional part of the number after taking the square root.

For example:

if (distanceSquared < radius * radius)
{
  (fill solid pixel)
}
else
{
   distance = sqrt (distanceSquared);

   if (distance < radius + 1)
   {
     alpha = 1 - (distance - radius);

     (blend pixel at alpha)
   }
}

This techinque is used in my open source library:

https://github.com/vinniefalco/LayerEffects.git


In addition to hotpaw2's answer:

  Presumably, the signal originated from a binary output on the transmitter end.

I'd say that is most probably wrong.

The typical digital transceiver based on complex baseband looks like this:

Binary data -(0/1)->
 Symbol mapper =(complex-valued symbols)=>
  pulse shaper =(complex-valued, "smooth" digital baseband signal)=>
   synchronous dual-DAC =(I & Q signals)=>
    filters =>
     I: mix with cos(2 pi f_RF) --\
                                   subtract--> filters -(RF signal)-> Antenna
     Q: mix with sin(2 pi f_RF) --/

On the transmitter end, you'll see some RF signal – and that is always analog, never binary.

  I haven't found any helpful documentation on how to do the final piece of this.

You cannot have read very much, because every textbook I know introduces both, the direct conversion transmitter and the direct conversion receiver (note: these go by different names, also zero-if or quadrature receivers).

Since you're already getting complex samples, we may assume you're working with a receiver that works like this

Antenna -(RF signal)->
 band-pass filter ->
  multiply with cos(2 pi f_RF) -> low-pass filter -(I signal) ->
  multiply with sin(2 pi f_RF) -> low-pass filter -(Q signal) ->
   =(IQ signal pair)=> synchronous dual ADC =(digital complex baseband signal)=>
    PC software

So, the analog part has been done for you.

Now, it's up to your PC software to do (not necessarily in this order, but all of these are typically necessary):

estimation/correction of the difference in f_RF (no two physical oscillators are perfectly identical)
estimation/correction of the difference in phase (i.e. the cos() of receiver and transmitter didn't "start" at the same time, and the signal was "in the air" for some non-zero amount of time due to limited speed of light)
reversing the pulse shaping of the transmitter
estimation/correction of the symbol rate (ie. how many digital samples is a single symbol long?)
guessing which symbols got sent by the transmitter
converting these symbols back to the bits they stood for

In reality, there's a few additional steps (channel coding, hardware imperfection corrections, equalizing, …).

Anyway, the order of these steps, and how exactly they're done depends completely on what kind of modulation scheme/waveform the transmitter produces; notice how you can't use a 3G mobile phone to get the bits your cordless landline phone sends, or use your WiFi card to receive digital TV, or use your CCTV receiver to receive satellite telephony.

Even for receivers for the same kind of signal, there's a lot of choices a SDR engineer makes – the type of timing recovery mechanism, for example, acceptable computational complexity when converting samples to symbols vs error rate when doing that, and a lot, lot more.

Thus,

  I realize this is most likely extremely simple,

Nope. It's the hardest problem, and pretty much the main problem, of  digital communications technology, the science backing one of the biggest technical sectors of modern days. Building good receivers is hard, in general. There's a lot of tools to help you deal with that. But: even these tools require solid understanding of the theory behind digital communications before you can go out and write your own SDR receiver.

But don't despair:

There's a lot of "simple" systems that are relatively easy to receive, and a tutorial alone might get you pretty far. It's an exciting field of technology, and I wouldn't have spent a solid  part of my life on that if it wasn't fun!

I've witnessed literally dozens of times how being able to actually get binary data out of something that is just a digital representation of an analog signal makes people a) happy and b) makes the understanding of the theory take giant leaps.

I'll encourage you to experiment with digital signals!

Based on your assumption, the transmitted signal was binary, I'll go ahead and take an audacious guess at what your profession is: Is it possible you are a Computer Science student?

I can see this signal world might be a bit daunting if you're coming from a background where you can safely assume that for

input---wire/file/network-->output

input==output, but that's simply not the case for digital communications.

In any case, no matter what you are, I'd recommend reading a bit on signal theory. The GNU Radio wiki has a page on recommended literature. Get one of the basic DSP books, and one of the basic Digital Comm books; some of them are freely available as PDF.

Parallely, take a 5-minute read of the introduction chapter of the GNU Radio guided Tutorials. These tutorials are meant to slowly, and carefully, and sequentially (don't skip a chapter) allow you to build a simple digital receiver. What they do is teach a tool, but they need to rely on some signal theoretic background, so they might be not that perfectly clear if you lack that.


If all you are looking to do is to take the 2D grayscale image and turn it into a 3D representation for visualisation purposes, then perhaps you could consider a function that accepts the image and produces a VRML (or X3D) string describing your 3D model. That model would then be passed to a standard java graphics toolkit (for example Java 3D) which would handle the rendering for you.

Here is a simplistic example:
For the 2x2 grayscale image of X = [0,1;1,0], your VRML indexedFace node would look like this:

IndexedFaceSet {
       coord Coordinate {
         point [-1 1 x(0,0), 1 1 x(0,1), 1 -1 x(1,1), -1 -1 x(1,0)]
       }
       coordIndex [ 0 1 2 3 ]   # That is, just one face
     }

In the above snippet, "x(i,j)" is a placeholder for the corresponding value from the image (which will be formated to the actual number) and as you can see it represents the Z dimension of your indexedFaceSet.

Actually, the structure of "point" is just a list (in row major order) of all of you image points, which is very easy to generate by iterating through all pixels of the image. Similarly, the structure of coordIndex is also set to quads (faces defined by 4 points in space) which map to the image's grid. The indexedFaceSet provides capabilities to define the colour of each vertex (each defining point of the face) as well in a way similar to the way "point" is defined.

More information about IndexedFaceSet can be found at:

VRML: http://www.cs.vu.nl/~eliens/documents/vrml/reference/CH3.HTM#IndexedFaceSet

X3D: http://www.web3d.org/documents/specifications/19775-1/V3.2/Part01/components/geometry3D.html#IndexedFaceSet

To render the model in Java, please see https://stackoverflow.com/questions/16625858/plug-in-display-3d-objects-from-vrml-or-x3d

Regarding the detection of spots, if these are images from electrophoresis, please see chapter 3 in this thesis for a comprehensive set of ways of detecting the blobs and extracting their main features: http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.86.4325&rep=rep1&type=pdf

I hope this helps


Lagrange Multipliers Optimization - Complex Functions

I found the following in Charles Therrien's "Discrete Random Signals and Statistical Signal Processing" in one of the Appendicies.

Say you have the function $Q(a)$ you wish to minimize such that $C(a)=0$, where $C(a)$ may be complex valued and $a$ may be a complex vector. The constraint really represents two real-valued constraints.
$$C_r(a)=0,\qquad
C_i(a)=0,
$$
where $C_r(a)$ and $C_i(a)$ represent the real and imaginary components. The Lagrangian can be written as:
$$\mathcal{L}=Q(a) +\mu_1C_r(a)+\mu_2C_i(a),
$$
where $\mu_1$ and $\mu_2$ are real valued.

This can also be written more compactly using a single complex parameter $\lambda=\lambda_r+j\lambda_i$ as
$$\mathcal{L}=Q(a) +\lambda C(a)+\lambda ^*C^*(a)
$$
where $\lambda^*$  and $C^*(a)$ denote the complex conjugate of the quantities. The reason for this is that
$$
\lambda C(a)+\lambda ^*C^*(a)=2Re[\lambda C(a)]=2\lambda_rC_r(a)-2\lambda_iC_i(a)
$$
and letting $\mu_1=2\lambda_r$ and $\mu_2=-2\lambda_i$, then the two Lagrangians  are seen to be equivalent.

Note - to apply the gradient to a real valued function with respect to a complex parameter you should also be familiar with Brandwood "A complex gradient operator and its application in adaptive array theory" IEE Proceedings 130(1) 11-16, Feb, 1983. This essentially treats the variables $a$ and $a^*$ as two separate variables and you only have to take the derivative wrt to one of them to perform a minimization.


phase offset extraction using matched filter

This is a common implementation for carrier recovery and tracking.

A matched performs a correlation of the received sequence with the expected sequence as given by:

$$\sum_{n=0}^{N-1}r[n]x^*[n]$$

Which states a correlation over N samples is performed by multiplying and accumulating the time-aligned received sequence $r[n]$ with the (complex conjugate of for complex signals) the expected ideal sequence $x[n]$. Note how any phase shift on $r[n]$ would result in the same phase shift on the result (any phase introduced on the ADC sampling clock would be translated to $r[n]$ directly as the sampled received signal):

$$\sum_{n=0}^{N-1}r[n]e^{j\theta}x^*[n] = e^{j\theta}\sum_{n=0}^{N-1}r[n]x^*[n]$$

Often a frequency offset is involved (which is a phase linearly changing with time), in which case we need to be aware of the resolution bandwidth of the matched filter, which is how the magnitude of the correlation result reduces as a function of frequency offset. For example If the matched filter is on constant magnitude coefficients (such as a moving average filter), the correlation magnitude will go down as a Sinc function with the first null at $1/T$ Hz where $T$ is the time length in seconds of the correlation.

But given the correlation result provides the carrier phase as shown above, if there is a small frequency offset, we can measure that precisely with successive results of the matched filter output given we will be observing a change of phase with a change in time (which is frequency given by $d\phi/dt$).

This is demonstrated below with captured data for live GPS signals where the matched filter (in this case for satellite SV24) is providing a sliding correlation over the 1 ms long GPS course-acquisition code, resulting in the spikes every 1 ms in the plots below indicating when the code aligns with the received signal. However the correlation is going up and down from to the underlying phase rotation due to a frequency offset.

This rotation is clearer on an IQ plot of the same result above, plotting the magnitude and phase on a complex polar plot:

By evaluating the angle over time of the correlation peaks, the frequency offset was determined in this case to be 38 Hz, and after correction the same above plots appear as follows, revealing the GPS data message (we see from the IQ plot that a small residual phase error exists that could also be removed):


Image processing: Image recognition and classification on fish

How can I get the specific key points?

I'm not sure but I think you are talking about Feature matching using SIFT.

Extract keypoints using SIFT from reference images of what you are trying to match.
Extract keypoints using SIFT from the test images.
Try to match the keypoints from reference images to test images using algorithms like RANSAC to get the position of the features(in your case eye,dorsal and pelvic fin).

You can find more information here.


What is the best way to detect rectangular pulse for AFSK demodulation?

How about you take a short time Fourier transform with window length equal to the granularity with which you want to locate the zero crossing. At a zero crossing the signal will have high frequency components. You could detect those by running a short time Fourier transform on the incoming samples to identify the location of zeros. For ex: if you window 5ms of data and take FFT and see if there are high frequency components of significant magnitude then this would mean that this window contains the zero crossing. If you want, you can go even narrower like 1ms, but you be would loosing on frequency resolution. So, pick an optimum resolution in frqeuency and required zero crossing accuracy in time and window the incoming signal and keep taking the FFTs. You could go for rectangular windows.

Even if noise is high it will be constant in spectrum (atleast close to it as is theoretically for white nosie) and this should not be a big trouble in FFT domain, unless it is higher compared to signal power itslef.

For just 13 sample long pulse durations use the following emperical algorithm

Maintain two numerical derivate with a reference point x samples apart. Fir ex: consider "sample number" $N_o$, take two derivates between $N_o$, $N_o -x$ call this gradient 1, and second derivative between $N_o$, $N_o + x$, call this gradient 2, take the magnitude of these two gradients, if they are greater than a threshold, usually gradient is close to 90 degree near the  zero crossing, then the point $N_o$ is a good approximation of zero crossing.  $x$, and the threshold can be tuned based on pulse length.


I cannot run your code because I have not matlab on this pc, but I can try to give you some advice.

The first thing I will check is the way you defined the cut-off frequencies. The Matlab function uses normalized frequencies (look at the examples here http://it.mathworks.com/help/signal/ref/buttord.html), which means you have to define it in this way: lets assume you want to cut your signal at 5 Hz and your sampling frequency is $Fs$

cut_freq = 5;
Wn = cut_freq/(Fs/2);

Of course if you have a band-pass filter you have to define 2 cut-off frequencies and therefore $Wn$ will be a vector of 2 elements.

Hope this helps. As soon as I get my laptop I will try to run your code and see what happens.


Pitch fluctuates unexpectedly in simple sine function

You are using 32-bit floats, and not resetting the phase by subtracting 2 pi.  That means the phase will eventually climb up to where the sin() function's phase unwrapping algorithm doesn't have enough bits of quantization, or valid mantissa, left over afterward unwrapping.

Add this:

if (phase > M_PI) { phase -= 2.0 * M_PI; }

inside your loop somewhere.


Detection of number of yarns in the fabric image

A fourier transform will give you the period of any periodic features in the image - eg it will tell you there is structure event 3pixels horizontally and 5 pixels vertically, from this and the pixel scale you can calculate the number of yarns


When to use Fourier, Laplace and Z transforms?

It's natural consequence of applying a transform to a convolution relation. The output $y(t)$ of an (continuous-time) LTI system is described by a convolution integral :

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

And when you apply a Fourier transform on this relation, it turns out to be a multiplication in the transform domain such as:

$$ Y(j\omega) = \mathscr{F}\{ y(t) \} = \mathscr{F}\{ h(t) \star x(t) \} $$
$$ Y(j\omega) = \mathscr{F}\{ h(t) \} \cdot \mathscr{F}\{ x(t) \}  = H(j\omega)\cdot X(j\omega) $$

Similary you can also apply a Laplace transform on it :
$$ Y(s) = \mathscr{L}\{ y(t) \} = \mathscr{L}\{ h(t) \star x(t) \} $$
$$ Y(s) = \mathscr{L}\{ h(t) \} \cdot \mathscr{L}\{ x(t) \}  = H(s)\cdot X(s) $$

And exactly the same happens for discrete-time LTI systems.


Basic Fourier Series Not Understanding how to interpret n values

I'm assuming the "n" you are referring to is the coefficient index, also referred to as the harmonic number; i.e. the n=3 term refers to the "third harmonic".

The n=0 term is the average value of the function, and is sometimes called the "DC component" because it does not depend on time.

Depending on the version of the Fourier series you are using, the terms in the series could go from negative infinity to positive infinity, or they could just go from n = 0 to positive infinity. The trigonometric version of the Fourier series (with sine and cosine) starts at n=0 and goes to infinity, whereas the complex exponential version (i.e. $c_n e^{j\omega_0 n t}$) has terms for negative and positive values of n.

The two versions are equivalent because the negative n terms cancel out the imaginary parts of the positive n terms in the complex exponential series.


Checking for VHF pulse with SDR and python

As far as the big picture goes, I assume you want to do both optimal detection of pulses and then tracking of the transmitters based on the detections you receive.  This answer will only deal with taking steps to optimize detection using a low cost RTL-SDR, and then I'll address some of your specific questions in your post.

In AWGN, the optimal detection filter is a matched filter, also called a correlator or a pulse compressor.  Creating a matched filter is conceptually very easy.  If $p[n]$ is a sequence of $N$ samples representing the pulse you are trying to detect, the matched filter $h[n]$ is the time-reversed, complex conjugate of the pulse, so $h[n] = p^*[(N-1)-n]$.

Let's take a look at how your SDR is going to present received pulses to you (assuming, for now, you haven't done any downsampling):

$F_s = 2.6 \space\mathrm{Msps} \\
t_p = 10 \space \mathrm{ms}\\
f_p = 50 \space \mathrm{kHz}\\
N = t_p \cdot F_s = 26,000 \space \mathrm{samples}$

$$p[n] = A e^{j\left(2\pi \frac{f_p}{F_s} n + \phi_0\right)} \quad \mathrm{for} \space n \in \{0, 1, \ldots, N-1\}$$

Thus the matched filter would be

$$h[n] = A e^{-j\left(2\pi \frac{f_p}{F_s}(N-1-n) + \phi_0\right)} \quad \mathrm{for} \space n \in \{0, 1, \ldots, N-1\}$$

Ideally, when filtering the ideal incoming pulse $p[n]$, when the pulse and the filter exactly line up, the filter output peaks at

$$\max_n y[n] = \max_n \left(p[n]*h[n]\right) = A^2\sum_{k=0}^{N-1} |e^0|^2 = A^2N$$

That peak should be well above the original received pulse and the filtered AWGN.

Now for the first bits of bad news.  The RTL-SDR crystal is not very precise or stable, so you won't precisely know $f_p$ ahead of time, nor will $F_s$ be exactly what you specify in real life.  The transmitter on the wildlife tag is also not very precise, which also contributes to $f_p$ being mismatched.  The transmitter on the wildlife tag is independent of your receiver, operating at some unknown distance and power-level, so you won't have $A$ as the received amplitude and the initial phase you receive will not be $\phi_0$.

The good news is that you can use a PLL to estimate and track $f_p$ and $\phi_0$, and also effectively correct for $F_s$ not being exact.

Normally one would use the estimates from the PLL to correct the incoming signal and then pass it to the matched filter to create peaks for detection.  However, since in your special case, your pulse is an unmodulated tone, you can take a shortcut and use the output of the PLL to implement an adaptive matched filter, instead of the normal steps of performing correction and then passing through a fixed matched filter.  This is what is going on in this section of a GNURadio flowgraph pictured below:

The "PLL Ref Out" block locks on to the strongest tone in the incoming signal and outputs a phase locked tone.  The "Multiply Conjugate" and "Moving Average" blocks implement the matched filtering operation, since the "Moving Average" block with a scale factor of 1 is just rolling summation.  The "Complex to Mag" block eliminates the effect of any residual phase offset from the "PLL Ref Out" block.

Below is a plot of the matched filtered output vs. the input data stream, showing the matched filter peaks from two pulses which are about 1.5 seconds apart:

The following picture is a zoom of the second pulse, to show the input pulse:

Once you have these matched filtered/correlation peaks, you can use whatever adaptive thresholding/outlier detection technique you want to detect the peaks from the noise.  I'll suggest Consecutive Mean Excision on the time domain data, if you don't have a favorite adaptive thresholding algorithm.  You can also use a high level "PLL" to track pulses that you think will be about 1 seconds apart to track a single transmitter.  The "PLL Ref Out" block I use also provides a frequency estimate output, so you can distinguish between transmitters at slightly different frequencies.

OK, now to step back an take a look at the whole GNURadio flowgraph I use:

Especially the filtering and downsampling stages:

To get the best noise performance from SDRs, it is best to use the highest sample rate the SDR provides, and then filter and downsample to the narrowest bandwidth possible that still contains your signal of interest plus any allowance for frequency offset errors.

Note that I have three stages of decimating filters and one final channel filter.  I broke the filtering and downsampling down into 3 stages of /4, /8, and /8 to keep FIR filter lengths and CPU utilization manageable.  Filtering and downsampling by /256 in one stage would have resulted in an extremely long filter and consumed much more CPU.  I have a fourth filter stage just to eliminate as much noise as I could, and still be sure to capture the pulse on one particular channel of the channels that the wildlife tags can use.  Note that the frequency translating FIR filter in the first stage can be used to shift to a different channel frequency within the band used by the wildlife tags.

So now for final bits of bad news.  SDR RF frontends are horrible in that they are wide open with very broad analog filters, if they have any at all.  They let in all sorts of RF interferers, some of them strong enough to clip inside the amplifiers and create intermodulation products that effectively raise your noise floor.  Cheap SDR units are not shielded from EMI and often bring in noise or interference from nearby electronics, such as the computer they are plugged into and its USB or Ethernet interface.  RTL-SDR units have gains in the RF frontend that can affect your noise performance -- RF Gain, Mixer Gain, IF Gain, and Baseband Gain -- and you are never provided enough information to set them properly for best performance.  The auto gain setting is usually suboptimal.  The LO leaks into the RF or IF front end and mixes in showing up as a mixing product at DC.  The direct downconversion can have I/Q imbalance, cause "ghost" signals reflected about DC from real signals in the spectrum.

If you want the best performance, you have to deal with the above problems with the RF frontend.  I recommend:

Get a higher quality SDR.  RTL-SDR units will always be poor performers.
Use a USB extension cable and get the RTL-SDR away from the computer and any nearby sources of VHF interference.
Use copper foil tape connected to a ground to shield the RTL-SDR unit, which is normally only in a plastic case.
Use a directional antenna.
Add an external LNA and SAW filter to limit the amount of analog RF bandwidth going in to the front end.
Set the mixer gain to the nominal mid level and leave it alone as long as the unit appears to be downconverting signals properly.
Set the RF Gain as high as you can, before clipping your signals of interest in the IF stage.
Set the IF Gain stage to its lowest setting at first and then increase until just before clipping your signals of interest in the baseband stage.
Don't tune the RTL-SDR exactly to your frequency of interest, but offset by at least a few kHz, to avoid the DC spike corrupting your signal of interest and pulling off your PLL.
Most RTL-SDRs that I've used can't support 2.6 Msps reliably.  I believe 2.4 Msps is the safe, maximum rate across all RTL-SDR variants.

I'll stop this rambling answer here.  I hope this helps.


Fast fourier transform and frequency domain spike in wrong area

You see one at index 8 because MATLAB's indexing is one-based, not zero-based as the equations that describe the DFT are. You get two spikes because your input signal is real and is therefore conjugate-symmetric in the frequency domain. If you work out the math, you'll see that the Fourier transform of a $\sin$ has two spikes symmetric about zero frequency. In order to get a single spike, you need to use a complex exponential function at the same frequency instead.


You are looking at an FIR filter whose transfer function is

$$
\begin{align}
H(z)&=(1-e^{j\omega_0}z^{-1})(1-e^{-j\omega_0}z^{-1})\\[10pt]
&=1-(e^{-j\omega_0}+e^{j\omega_0})z^{-1}+z^{-2}\\[10pt]
&=1-2\cos(\omega_0)z^{-1}+z^{-2}
\end{align}$$

So as you can also see in the comment, there are two zeros at $e^{\pm j\omega_0}$ that null the frequencies $\pm\omega_0$.

Obviously, here we have $\omega_0=0.44\pi$. The filter's response is

The impulse response is $$h[n]=\delta[n]-2\cos(0.44\pi)\delta[n-1]+\delta[n-2]$$

Assume an arbitrary input $x[n]$. The output $y[n]$ is given by:

$$y[n]=x[n]*h[n]$$ Therefore $$y[n]=x[n]-2\cos(0.44\pi)x[n-1]+x[n-2]$$


No it does not satisfy the condition. Simply take an example:

$$n = 1 \implies y[1] = x[2]$$

Hence the output value at the present time $n=1$ depends on a future value of the input at time $n=2$. This violates the causality principle.


FIR filter design and implementation meaning

I don't fully understand how you are expressing your question in English, so please accept my attempt.

Discrete convolution is
$$ y[n]=\sum_{k=-\infty}^{\infty} x[n-k] h[k]
$$
where $y[n]$ is the output at time index $n$ and $h[k]$ is the impulse response.  The index $k$ ranges over $-\infty$ to $\infty$ if $h[n]$ is an infinite impulse response.

A causal convolution, i.e. the filter can only have output $y[n]$ for $x[n]$ ,$x[n-1]$ ,$x[n-2]$ $\dots$ $x[n-\infty]$
$$ y[n]=\sum_{k=0}^{\infty} x[n-k] h[k]
$$
if $h[0]\ne 0$, $h[1]\ne 0$, $h[2] \ne 0$, $\dots$ $h[N-1]\ne 0$ , $h[N]=0$,  $h[N+1]=0$, $\dots$,  $h[\infty]=0$, is of finite length, it is a FIR filter and

  what is the output of the convolution meaning which is used in the FIR
  filter ?

$$ y[n]=\sum_{k=0}^{N-1} x[n-k] h[k]
$$

  how can conclude the output of the FIR filter in the output according
  to the input x[n]?

$y[n]$ depends on a finite sum of $x[n]$  .  A FIR filter also need not necessarily be causal, but it doesn't change that $y[n]$ only depends on a finite sum of $x[n]$

I don't understand how design has anything to do with your question.


Can I square signal in the frequency domain to reduce noise?

Squaring your signal is multiplying it with itself.

Multiplication in frequency domain is convolution in time domain.

So you're correlating the signal in time domain with its time-inverse. I don't think this has beneficial effects aside from a very specific class of signals (namely, the signals which are their own time-inverse).

You forget why you're interested in SNR in the first place. It's not for the sake of SNR itself. It comes from the assumption that your signal is essentially unharmed but overlaid with additive noise. You're breaking that assumption.

Also, the effect is purely optical; you know, squaring positive numbers is a monotonous operation. So if you need to find something "clearly larger" than something else, you can do this with the original graph just as well.

Adapting the Y-axis to be able to spot signals of very different magnitudes is a very common thing to do, however: for PSD-style plots you typically don't look at abs(spectrum) in a linear scale, but convert it to dB, which is a logarithmic scale; in fact, it hence does quite the opposite of your polynomial scale, but it's much more useful for real-world signals.

There's a few cases where you'd square the time signal (e.g. to allow for easy timing recovery of BPSK signals), but from the top of my head, I can't think of anything where you'd work with the squared DFT; might be useful if we know that a signal has a time-inverse part and were looking for that.

Note that Laurent correctly pointed out that "squaring" is a bit of an imprecise term, and in the complex sense you'd normally use it as "multiplying with the complex conjugate", but that doesn't make much of a difference here.


Subtract a digital filter from another filter

If (and only if) the filters (either FIR or IIR) are linear, then you can subtract their coefficients term-wise, supposing they are correctly aligned (below, on the $0$ index), treating them as infinite series with trailing zeroes when the coefficients are not defined. For instance with:

$$h = [h_0, h_1, h_2]$$
and
$$g = [g_{-1},g_0, g_1, g_2,g_3]$$

since the convolution with $h$ yields the same result as the convolution with a zero-padded $h$ (and $g$)

$$h' = [\ldots,0,0,h_0, h_1, h_2,0,0,\ldots]$$
the difference filter $d$ will have (symbolically $h'-g'$):

$$d = [\ldots,0,-g_{-1},h_0-g_{0}, h_1-g_{1}, h_2-g_{2},-g_{3},0,\ldots]$$

For FIR filters (or truncated IIR), you can limit to a finite sequence bounded by the smallest and the biggest index of both filters.


How can I detect potholes using an accelerometer mounted in a vehicle?

You are not really interested in the frequency information in your signal so much as detectable features - as such I would suggest that there isn't that much utility in an FFT here.

What you really want to do is correlate the continuous 'car signal' recorded in the car with a different signal which is the 'signature' of the pothole bump. You will need to acquire this signature by recording it in a car going over a pothole (perhaps recording many and using them to produce a generic signal which is a good representation of a pothole). Once you have it you can cross correlate it with your car signal. You may want to do this in the frequency domain in which case you should cross correlate the FFT of your signature with a sliding FFT of the car signal.

You may want to process the car signal before you use it to remove, for example, high frequency noise (with a lowpass filter). You can do this pretty effectively by eye - play around with various filters on car signals to see which removes the most extraneous information while leaving the important pothole information intact.


PCA to reduce dimensionality to 99% variance

If what you want is PC compression, a more useful form of PCA would be achieved via Singular Value Decomposition (which is, in many cases, more accurate and faster than an eigendecomposition).

Assuming that X is a "tall" matrix with dimensions [n, p], n>p;
Then:

 [U, S, V] = svd(X);

eig_vals = diag(S.^2); %recall that eigenvalues of cov(X) are the singular values of X, squared.
percentages = cumsum(eig_vals)./sum(eig_vals); % cumulative variance.

% Find the 99% variance
PC99 = find(percentages >= 0.99, 1);

X_reduced = U(:,PC99)*S(PC99,PC99)*V(PC99,PC99)'; % compressed version of X

If X is a wide matrix [n < p], you have to adjust the above expression a bit (all of V will be used, but only part of U, unlike here where V is truncated as well). We leave it as an exercise to the interested reader!

I will leave wiser people to comment on the judiciousness of PCA for compression. In many cases there are much more efficient schemes, such as wavelet decomposition followed by "thresholding" (zeroing out the coefficients  that fall below a certain threshold). But since your question assumes that you know what you want, so does the answer.


why we do multiple aliasing of high order harmonics

Aliasing is more like circular wrapping around rather than just mirroring.  Frequencies wrap circularly between -Fs/2 to Fs/2, an 800 Hz range from -400 to 400 at a sample rate of 800 sps.
120 is inside -Fs/2 to Fs/2, so doesn't need to wrap around at all.
But 480 wraps 80 past 400, and shows up at -320, which looks like 320 if you don't care about phase.
And 1120 wraps all the way around the circles range of 800 Hz, plus another 320, so it shows up at 320.
12000 Hz (12 kHz) would wrap around 15 times and show up identical to DC or 0 Hz.


Inverse Continuous Wavelet Transform derivation?

Summary: the dual wavelet's role is analogous to that of $e^{j\omega t}$; it undoes the wavelet's convolving with the signal (integrated inner product). The main intricacy's indeed in normalization; there's no exact answer without digging in math (done below) - but intuitively, the seemingly matter-worsening $a^{-2}$ is due partly to the dilating nature of the wavelet, i.e. $\psi(t/a)$, which upon inversion requires downscaling not only to properly preserve energy, but reconstruct correctly by undoing the inner product.
Below are more "proofs" rather than "derivations", as they work backwards to show equality with $f(t)$, but I term them latter as full-fledged proofs are beyond the scope of this answer.

Derivation 1: Convolution Theorem
Pasted for convenience the CWT and iCWT
$$
W_f(a, b) =   \int_{-\infty}^{\infty} f(t) \frac{1}{|a|^{1/2}} \psi \left(\frac{t - b}{a}\right) dt \tag{0}
$$
$$
f(t) = \frac{1}{C_{\psi}} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} W_f(a,b) \frac{1}{|a|^{1/2}} \psi \left(\frac{t - b}{a}\right) db \frac{da}{a^2}, \tag{1}
$$
Begin by rewriting the CWT as a convolution: $W_f(a, b) = f \star {\bar\psi_a}(b)$, where $\psi_a(t)=a^{-1/2}\psi(t/a)$. Inserting into $(1)$, and letting the right side equal $b(t)$,
$$
\begin{align}
b(t) &= \frac{1}{C_\psi} \int_0^{\infty} W_f(.,a) \star \psi_a(t) \frac{da}{a^2} \\
     &= \frac{1}{C_\psi} \int_0^{\infty} f \star {\bar\psi_a}(b) \star \psi_a(t) \frac{da}{a^2} \tag{2}
\end{align}
$$
where we also expressed the $db$ integral as a convolution, and $"."$ indicates the variable over which convolution is computed. We prove $b=f$ by showing their Fourier Transforms are equal. With $\hat{} = \mathcal{F}$,
$$
\begin{align}
{\hat b}(\omega) &= \frac{1}{C_\psi} \int_0^{\infty}
    {\hat f}(\omega) \sqrt{a} {\bar{\hat \psi}}(a\omega) \sqrt{a} {\hat\psi}(a\omega) \frac{da}{a^2} \tag{3a} \\
                 &= \frac{{\hat f}(\omega)}{C_\psi} \int_0^{\infty} \tag{3b}
    |{\hat \psi}(a\omega)|^2 \frac{da}{a}
\end{align}
$$
With change of variables $\xi = a\omega$, we thus 'prove'
$$
{\hat b}(\omega) = {\hat f}(\omega) \left(\frac{1}{C_\psi} \int_0^\infty \frac{|{\hat \psi}(\xi)|^2}{\xi} d\xi \right) = {\hat f}(\omega). \tag{4}
$$
Interpretation:
It's tempting to attribute the $a^{-2}$ to the wavelet dilations, $\psi(t/a)$, as they manifest as rescalings of Frequency-domain components ($\cdot \sqrt{a}$), and while that's true, it's not the crux.
Recall when I said it's not just about energy? Consider what happens if the normalization is anything other than $a^{-2}$. To arrive at $(4)$, we relied on $a^1$ in denominator in $(3b)$, else there is no cancellation with $C_\psi$, so no reconstruction. Thus, one can conclude that $a^{-2}$ is set in retrospect - i.e. "what do we make it to cancel $C_\psi$".
Note that we can't go the other way, i.e. redefine $C_\psi$ such that $a^{-2}$ is no longer needed - at least not entirely; per above derivation, it's clear that $\psi$ must somehow drop, thus $C_\psi$ must include it. The usage of exactly $C_\psi$ might be convention-based, since $C_\psi$'s applicability spans much outside the iCWT as the fundamentally defining criterion of wavelets (Farge, 1992).
So how to invert L1-normed CWT? Obvious from above.
Notes:

Proof followed Mallat, Theorem 4.4
Split integrals based on identity
Assumed $f \in \mathbf{L^2}(\mathbb{R})$ (square-integrable)
$\psi$ was also assumed real, and note the tilde is missing. I'm unsure why such a restriction was made, or whether Mallat implies that all real $\psi$ are self-dual, but the rescaling identities (e.g. $a \Rightarrow a^{-1}$ via $\mathcal{F}$) should hold regardless. Derivation 2 does not restrict $\psi$ to real.

Derivation 2: Resolution of Identity
For all $f, g \in \mathbf{L^2}(\mathbb{R})$ (real or complex), below holds:
$$
\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \frac{da db}{a^2} W_f(a, b) \overline{ W_g (a, b)} = C_\psi \left< f, g \right> \tag{5}
$$
Begin proof by transforming to frequency domain:
$$
\begin{align}
(5) = \int \int \frac{da db}{a^2} & \left[
      \int d\xi {\hat f}(\xi) |a|^{1/2} e^{-jb\xi} \overline{ {\hat \psi}(a\xi)} \right] \\
      \cdot & \left[ \int d\xi' \overline{ {\hat g}(\xi')} |a|^{1/2} e^{j b\xi'} {\hat \psi}(a\xi')
      \right] \tag{6}
\end{align}
$$
This is similar to derivation 1, except $\psi$ replaced with another CWT of some different function $g$. In the brackets we have the Fourier transform of $F_a(\xi) = |a|^{1/2} {\hat f(\xi)} \overline{\hat \psi (a\xi)}$ times the Fourier transform of complex conjugate of $G_a(\xi) = |a|^{1/2} {\hat g(\xi)} \overline{\hat \psi (a\xi)}$. By the unitarity of the Fourier transform, we have
$$
\begin{align}
(6) &= 2\pi \int \frac{da}{a^2} \int d\xi F_a(\xi) \overline{G_a(\xi)} \\
    &= 2\pi \int \frac{da}{|a|} \int d\xi {\hat f(\xi)} \overline{ {\hat g}(\xi)} |{\hat \psi}(a\xi)|^2 \\
    &= 2\pi \int d\xi {\hat f(\xi)} \overline{ {\hat g}(\xi)}  \int \frac{da}{|a|} |{\hat \psi}(a\xi)|^2 \ \ \ \ \text{(Fubini's theorem)}  \\
    &= C_\psi \left< f,g \right>,
\end{align}
$$
where in last step a similar change of variables was used, $\zeta = a\xi$. $(6)$ can be read as
$$
f = C_\psi^{-1} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \frac{da db}{a^2} W_f(a, b) \psi^{a, b}, \tag{7}
$$
quoting the source,

with convergence of the integral "in the weak sense", i.e., taking the inner product of both sides of $(7)$ with any $g \in \mathbf{L^2}(\mathbb{R})$, and commuting the inner product with the integral over $a, b$ in the right-hand side, leads to a true formula.

This seems to say that the original derivation was for the inner product $\left< f, g \right>$, whereas $(7)$ reads only for $(6)$, which will have somewhat different properties. The "original" is recovered by convolving with a $g$. -- I'm also unsure what exactly "in the weak sense" means; this is the most I found (clarifications welcome).
Once again, it is evident that $a^{-2}$ was set in retrospect, depending on how the forward transform $W$ was normalized.
Source: Daubechies, Ten Lectures on Wavelets, Chapter 2


Alas, optical flow is a difficult problem too ;-)

Well, to be more constructive, here are a few algorithms that should be worth trying (or have been tried on this particular sequence) :

re-train your bags of features on a databse of vehicles more representative (in size and orientation) to your actual problem in order to obtain better results
use the fact that the ground is a flat plane to do some parametric optical flow (search for an affine flow) or to compute some affine registration between the frames of the sequence. The moving vehicles will then be outlier from this dominant motion
use some optical flow algorithm to compute the flow, then try to classifiy / cluster the optical flow vectors (this is still a widely open problem!). Depending on the language you use, you can use OpenCV's optical flow, the one from TU Graz, D. Sun's optical flow, or even mine ;-). Note however that segmenting the flow will be a non-trivial task that you should probably do in two steps: global (dominant) motion etsimation, then small motion detection.


You've pointed out a very important distinction between theory and practice. In theory, as suggested by your book and in Fat32's answer, modulation schemes where all information resides in the phase of the signal, not in the amplitude, are called "constant envelope modulation". However, in practice our systems have finite bandwidth and instantaneous phase jumps cannot be realized. Consequently, in such cases the envelope will have dips where (theoretical) phase jumps occur.

Nevertheless, there are systems which really have a constant envelope. Such systems produce a (complex baseband) signal of the form

$$s(t)=Ae^{j\phi(t)}\tag{1}$$

where $\phi(t)$ is a continuous function of time. Such signals are well suited for extremely non-linear channels (amplifiers), and they are also more spectrally efficient than signals with a (theoretically) discontinuous phase.

Modulation schemes producing a signal as given by $(1)$ are called continuous phase modulation (CPM). Examples of CPM are continuous-phase FSK (CPFSK), and minimum-shift keying (MSK) and its variants (e.g., GMSK)