Would anyone be able to help out with the following discrete convolution question?

Hint

The simplest way is to use the Z transform property "convolution in time domain is multiplication in z domain".
See Z transform convolution

$$\mathrm{Z}(x[n]*h[n]) = \mathrm{Z}(x[n]) \times \mathrm{Z}(h[n])$$

Then you just need to do inverse Z transform.

Some typical Z transforms (including what is useful for you) can be found at Z transform pairs

Be careful with region of convergence (ROC).

It seems a homework?


Purpose of using polyphase filter bank (PFB)

One advantage of a polyphase filterbank approach is, as you guessed, that you can control the frequency response of each channel. When using a DFT alone, you have limited control over the frequency band covered by each bin (characterized by a Dirichlet kernel in the unwindowed case, or by the frequency response of whichever window function you select). This is sufficient for a lot of applications.

In some other applications, however, you might want tighter control over the per-channel frequency response. Say you wanted to construct a DFT-based spectrum analyzer with very particular specifications (e.g. -3 dB response at the midpoint between output bins, -80 dB response at a one-bin spacing from center). You can utilize the polyphase filterbank structure to implement whatever filter is needed to achieve the specified level of performance.

Another application is in cases where reconstruction is needed after the channelizer: if you're careful in the design of the polyphase filter, you can actually perform straightforward spectral modification (similar to the "ideal filtering" that many signal processing beginners attempt using the DFT) in the frequency domain, then use a synthesis filterbank to take the composite signal back into the time domain. This structure, with cascaded analysis and synthesis stages, is known as a transmultiplexer.


Relation between probability of a signal and its magnitude spectrum

Given a signal y(t), the magnitude spectrum  is
$$Y(f)=\left\lvert\int y(t)e^{-2\pi ft}\,\mathrm dt\right\rvert,$$
in other words the absolute value of the FT of the signal.

But that's only applicable for deterministic signals¹. If your signal itself is random, then you can't know $y(t)$ – you only know some stochastic properties of it. There's no $y(t)$ to transform if you don't actually know $y(t)$. For random signals, $y(t)$ can't be given, otherwise the signal wouldn't be random. You might observe a realization, but that realization doesn't describe the full stochastic process!
The definition of the power spectral density for random signals like your $y$ is the Fourier transform of the autocorrelation function. For that to be sensibly defined, you need your autocorrelation function to depend only on one variable, the time shift, which means your random signal needs to be wide-sense stationary. But, if that's the case, then the autocorrelation function (ACF) is defined as
$$\phi_{yy}(\tau) := \mathbb E \left\{ y(t) y^*(t+\tau) \right\},$$
and then the power spectral density, which gives your magnitude spectrum is its Fourier transform, $\Phi_{yy}(f) :=  \int_{-\infty}^{\infty} \phi_{yy}(\tau) e^{-i2\pi f\tau}\,\mathrm d\tau$.
Now, in the expectation $\mathbb E \left\{ y(t) y^*(t+\tau) \right\}$, there's your random signals probability density function:
\begin{align}
\phi_{yy}(\tau) &:=
\mathbb E \left\{ y(t) y^*(t+\tau) \right\} \\
&=
\int_{-\infty}^\infty
  \int_{-\infty}^\infty
    y'_1\,y'_2
    \,\cdot\,
    f_{
      (y({t_1}),y({t_1+\tau}))
    }
   (y'_1,y'_2)
    \,
  \mathrm d y'_2
  \;
\mathrm d y'_1
\end{align}
Notice that here, the density is a function of two values, and parameterizes over the time shift.

¹ you're also missing imaginary unit in the exponent and the integration bounds, but that's not the problem here. But, really, make sure to be able to write down the Fourier transform correctly.


Is the Dirac delta (impulse) signal a power signal or an energy signal?

[Added a reference on Schwartz's impossibility theorem for products of distribution]

The continuous Dirac delta  $\delta$ is not considered a true function or signal, but a distribution. From its wikipedia page:

  The delta function can also be defined in the sense of distributions
  exactly as above in the one-dimensional case.[25] However, despite
  widespread use in engineering contexts, (2) should be manipulated with
  care, since the product of distributions can only be defined under
  quite narrow circumstances.

It can be defined such that, for any function $f$ satisfying some important properties, and for $a\in \mathbb{R}$:

$$ \int f(t)\delta(t-a)dt=f(a).$$

From the best of my knowledge, those important properties are not satisfied by $\delta$, so one cannot directly replaced $f$ by $\delta$ and get a meaningful result. As far as know, the product of two Dirac distributions is not well defined, unless one talks about $n$-dimensional versions, or so-called "formal" manipulations as used in physics for instance, or more complicated maths.
A short account of Simplified production of Dirac delta function identities is provided by Nicholas Wheeler. If one wants to dig deeper, I'd suggest The Colombeau theory of generalized functions, by Ta Ngoc Tri, 2005:

  Soon after the introduction of his own theory, L. Schwartz published a
  paper in which he showed an impossibility result (see [Sch54]) about
  the product of two arbitrary distributions.

One result is  Schwartz impossibility result. It (somehow) says that if one wants to encompass the derivative of continuously differentiable functions whle keeping  Leibniz's rule of derivation, then one get $\delta^2(|x|)=0$.

However, from an informal point of view, sometimes used in DSP (and in physics), this "product" is, as far as I know, neither energy nor power.
From a logical point of view though, if it does not exist, one could affect this "product" a lot of properties...

Some related posts:

What is product of delta function with itself?
What is the square root of the Dirac Delta Function?
Products and compositions with the Dirac delta function, C. K. Raju, Journal of Physics A: Mathematical and General, Volume 15, Number 2


Have you had a look at the diphone synthesizer resources in the Festival Speech Synthesis framework?

The state of the art for Text2Speech (though still not practical for production due to GPU compute requirements) is Deepmind's WaveNet.

Their results are definitely smoother than concatenative or parametric methods.  Good luck!


Given 1 recording of my voice, how can I process it to sound like 2 independent recordings of my voice played in unison?

To make it sound natural you typically introduce some small random pitch, speed & gain modulations. That's not a trivial amount of work, especially if you want something that sounds natural and good and maintains the original phrasing.

This is a pretty common plug-in in audio processing. It's typically called "vocal doubler" or "voice doubler" and there are many off-the-shelf options available. Here is an example: https://www.izotope.com/en/products/vocal-doubler.html. No endorsement or commercial interest, but you can't beat the price :-)


What is the relation between eigenvalues and state-space response in control systems?

Just because eigenvalues apply in linear algebra doesn't mean that the same meaning applies everywhere. Think of the Fourier transform: does it only apply to the study of heat transfer? Same here. The A matrix is the companion matrix, which means its eigenvalues represent the roots of the system. For example:
$$\begin{align}
H(s)&=\dfrac{3}{s^2+2s+3}\tag{1} \\
\mathbf{A}&=\begin{bmatrix}{{0}\quad{-3} \\ {1}\quad{-2}}\end{bmatrix}\tag{2} \\
s_{1,2}&=-1\pm j\sqrt{2}\tag{3} \\
\end{align}$$
I'm afraid I don't quite understand what you mean by "from a non-mathematical (physical) point of view". The eigenvalues are the mathematical roots of the system, they tell about the stability, the time and frequency response. There's nothing non-physical about it. The system itself is physical -- as in it has a meaning -- even if only described on paper, such as the example above. It's just that you don't have to think about them from a single perspective, only.


Detecting direction of sound using several microphones

To extend Müller's answer,

  Should the microphones be placed in separate tubes in order to improve separation?

No, you are trying to identify the direction of the source, adding tubes will only make the sound bounce inside the tube which is definitely not wanted. The best course of action would be to make them face straight up, this way they will all receive similar sound and the only thing that is unique about them are their physical placements which will directly affect the phase. A 6 kHz sine wave has a wavelength of $\frac{\text{speed of sound}}{\text{sound frequency}}=\frac{343\text{ m/s}}{6\text{ kHz}}=5.71\text{ mm}$. So if you want to uniquely identify the phases of sine waves up to 6 kHz, which are the typical frequencies for human talking, then you should space the microphones at most 5.71 mm apart. Here is one item that has a diameter that is less than 5.71 mm. Don't forget to add a low pass filter with a cut-off frequency at around 6-10 kHz.

Edit

I felt that this #2 question looked fun so I decided to try to solve it on my own.

  Can phase correlation be calculated between 3 sources simultaneously somehow? (i.e. in order to speed up the computation)

If you know your linear algebra, then you can imagine that you have placed the microphones in a triangle where each microphone is 4 mm away from each other making each interior angles $60°$.

So let's assume they are in this configuration:

       C
      / \
     /   \
    /     \
   /       \
  /         \
 A - - - - - B

I will...

use the nomenclature $\overline{AB}$ which is a vector pointing from $A$ to $B$
call $A$ my origin
write all numbers in mm
use 3D math but end up with a 2D direction
set the vertical position of the microphones to their actual wave form. So these equations are based on a sound wave that looks something like this.
Calculate the cross product of these microphones based on their position and waveform, then ignore the height information from this cross product and use arctan to come up with the actual direction of the source.
call $a$ the output of the microphone at position $A$, call $b$ the output of the microphone at position $B$, call $c$ the output of the microphone at position $C$

So the following things are true:

$A=(0,0,a)$
$B=(4,0,b)$
$C=(2,\sqrt{4^2-2^2}=2\sqrt{3},c)$

This gives us:

$\overline{AB} = (4,0,a-b)$
$\overline{AC} = (2,2\sqrt{3},a-c)$

And the cross product is simply $\overline{AB}×\overline{AC}$

$$
\begin{align}
\overline{AB}×\overline{AC}&=
\begin{pmatrix}
4\\
0\\
a-b\\
  \end{pmatrix}
×
  \begin{pmatrix}
2\\
2\sqrt{3}\\
a-c\\
  \end{pmatrix}\\\\
&=\begin{pmatrix}
0\cdot(a-c)-(a-b)\cdot2\sqrt{3}\\
(a-b)\cdot2-4\cdot(a-c)\\
4\cdot2\sqrt{3}-0\cdot2\\
  \end{pmatrix}\\\\
&=\begin{pmatrix}
2\sqrt{3}(b-a)\\
-2a-2b-4c\\
8\sqrt{3}\\
  \end{pmatrix}
\end{align}
$$

The Z information, $8\sqrt{3}$ is just junk, zero interest to us. As the input signals are changing, the cross vector will swing back and forth towards the source. So half of the time it will point straight to the source (ignoring reflections and other parasitics). And the other half of the time it will point 180 degrees away from the source.

What I'm talking about is the $\arctan(\frac{-2a-2b-4c}{2\sqrt{3}(b-a)})$ which can be simplified to $\arctan(\frac{a+b+2c}{\sqrt{3}(a-b)})$, and then turn the radians into degrees.

So what you end up with is the following equation:

$$\arctan\Biggl(\frac{a+b+2c}{\sqrt{3}(a-b)}\Biggr)\frac{180}{\pi}$$

But half the time the information is literally 100% wrong, so how.. should one.... make it right 100% of the time?

Well if $a$ is leading $b$, then the source can't be closer to B.

In other words, just make something simple like this:

source_direction=atan2(a+b+2c,\sqrt{3}*(a-b))*180/pi;
if(a>b){
   if(b>c){//a>b>c
     possible_center_direction=240; //A is closest, then B, last C
   }else if(a>c){//a>c>b
     possible_center_direction=180; //A is closest, then C last B
   }else{//c>a>b
     possible_center_direction=120; //C is closest, then A last B
   }
}else{
   if(c>b){//c>b>a
     possible_center_direction=60; //C is closest, then B, last A
   }else if(a>c){//b>a>c
     possible_center_direction=300; //B is closest, then A, last C
   }else{//b>c>a
     possible_center_direction=0; //B is closest, then C, last A
   }
}

//if the source is out of bounds, then rotate it by 180 degrees.
if((possible_center_direction+60)<source_direction){
  if(source_direction<(possible_center_direction-60)){
    source_direction=(source_direction+180)%360;
  }
}

And perhaps you only want to react if the sound source is coming from a specific vertical angle, if people talk above the microphones => 0 phase change => do nothing. People talk horizontally next to it => some phase change => react.

$$
\begin{align}
|P| &= \sqrt{P_x^2+P_y^2}\\
&= \sqrt{3(a-b)^2+(a+b+2c)^2}\\
\end{align}
$$

So you might want to set that threshold to something low, like 0.1 or 0.01. I'm not entirely sure, depends on the volume and frequency and parasitics, test it yourself.

Another reason for when to use the absolute value equation is for zero crossings, there might be a slight moment for when the direction will point in the wrong direction. Though it will only be for 1% of the time, if even that. So you might want to attach a first order LP filter to the direction.

true_true_direction = true_true_direction*0.9+source_direction*0.1;

And if you want to react to a specific volume, then just sum the 3 microphones together and compare that to some trigger value. The mean value of the microphones would be their sum divided by 3, but you don't need to divide by 3 if you increase the trigger value by a factor 3.

I'm having issues with marking the code as C/C#/C++ or JS or any other, so sadly the code will be black on white, against my wishes. Oh well, good luck on your venture. Sounds fun.

Also there is a 50/50 chance that the direction will be 180 away from the source 99% of the time. I'm a master at making such mistakes. A correction for this though would be to just invert the if statements for when 180 degrees should be added.


Output of IIR filter not matching the frequency response - Proximity to fs/4

From the octave documentation on cheby1.m you can see that the cut-off frequency in radians is given by pi*Wc, where Wc is the cut-off frequency input argument of cheby1.m. So for a function call [b,a] = cheby1(n,Rp,Wc), you should test the filter's gain at the cut-off frequency with an input sequence cos(Wc*pi*n). Furthermore, choose a sufficiently long input signal so that you can observe the filter's steady-state behavior.
Apart from that it's a matter of sampling the time domain function. Note that the samples will generally not occur at the maxima of the sinusoid. I.e., the maximum magnitude of the output signal
$$y[n]=|H(\omega_c)|\cos(n\omega_c+\phi(\omega_c))$$
will not necessarily attain its theoretical maximum value $|H(\omega_c)|$.


Source Localization with the BLUE Estimator

One problem with your code is that the range calculations should be

R(i) = sqrt((xs-x(i))^2$\color{red}{\bf +}$(ys-y(i))^2);

Rn(i) = sqrt((xn-x(i))^2$\color{red}{\bf +}$(yn-y(i))^2);

But that doesn't appear to help.

The thing that seems to be the issue is that your noise variance is way too high.  If I set it to var = 0.00000001 then I get good estimates.

The reason is because you're using this variance (or a standard deviation of 0.0001) to add noise to your H matrix. Looking at the values there:

you can see you're adding something with a variance of 1 (standard deviation of 1), to something of the order of $10^{-5}$.  That means your SNR (signal-to-noise ratio) is about -100dB.  I know of no algorithm that will help you at SNRs that low.

You left out that Kay defines
$$
R_i(x_s,y_s) = \sqrt{(x_s - x_i)^2 + (y_s - y_i)^2}
$$
and
$$
R_{n_i}(x_n,y_n) = \sqrt{(x_n - x_i)^2 + (y_n - y_i)^2}
$$
where $x_s = x_n + \delta x_s$ and $y_s = y_n + \delta y_s$.

Then the aim is to write $R_i$ in terms of $R_{n_i}$.

Wikipedia has an example of Taylor series in two variables, and the first order approximation is just:

$$
f(x,y) \approx f(a,b) + (x-a) \frac{df}{dx} (a,b) + (y-b) \frac{df}{dy} (a,b) + \ldots
$$

Then all you need to do is apply Wikipedia's expression to Kay's to get

$$
R_i \approx R_{n_i} + \frac{x_n-x_i}{R_{n_i}} \delta x_s + \frac{y_n - y_i}{R_{n_i}} \delta y_s.
$$


Calculating normalized edge frequencies for filter design in Matlab

This is from the mathworks documentation on cheby2

For digital filters, the stopband edge frequencies must lie between 0 and 1, where 1 corresponds to the Nyquist rate—half the sample rate or $\pi$ rad/sample.

For discrete-time signals, we use the normalized frequency $\omega$ in radians (per sample), defined as
$$\omega=\frac{2\pi f}{f_s}\tag{1}$$
where $f$ is the actual frequency in Hertz, and $f_s$ is the sampling frequency.
Matlab uses $(1)$ normalized by $\pi$, i.e., edge frequencies etc. are defined by W$=2f/f_s$. The value W$=1$ corresponds to $f=f_s/2$, which is the Nyquist frequency.
Example: The cut-off frequency is $f_c=500$ Hz and the sampling frequency is $f_s=2000$ Hz (samples/second). According to $(1)$, the corresponding normalized frequency in radians (per sample) is
$$\omega_c=\frac{2\pi f_c}{f_s}=\frac{\pi}{2}\tag{2}$$
Consequently, when calling a Matlab routine you would use Wc$=\omega_c/\pi=\frac12$.


Difference between rceps and cceps in MATLAB

RTFM? rceps says:

y = real(ifft(log(abs(fft(x)))));

while cecps says:

h = fft(x);
logh = log(abs(h)) + sqrt(-1)*rcunwrap(angle(h));
y = real(ifft(logh));


Is there a closed form for the average energy of a QAM constellation?

From $\bar{E}=\frac{(M-1)d^2}{6}$ as the average energy in $M$-QAM, $M$ a perfect square and a square constellation, we have that the total energy is
$$E_{\operatorname{tot}} = M\frac{(M-1)d^2}{6}$$
Now, since the $M$-QAM signal is the sum of two orthogonal $\sqrt{M}$-PAM
signals of equal energy, then, writing $N$ instead of $\sqrt{M}$,
we have each of the orthogonal signals carrying a total energy of
$$\displaystyle \frac{1}{2}\frac{M(M-1)d^2}{6}
= \frac{N^2(N^2-1)d^2}{12}.$$
Since there are $M = N^2$ constellation points, the average energy of
a $N$-PAM signal with spacing $d$ is
$\displaystyle \frac{(N^2-1)d^2}{12}$.

Now, if $b$ is odd, a $2^b$-QAM signal constellation is
a rectangular grid with $2^{(b-1)/2}$ points in one direction and
$2^{(b+1)/2}$ points in the other.  Thus, the $2^b$ points contribute an
average of
$$\frac{(2^{b-1}-1)d^2}{12} +  \frac{(2^{b+1}-1)d^2}{12}
= \frac{(2.5M-2)d^2}{12}$$

Note that it is possible to use constellations other than rectangular when
$b$ is odd. For example, the $32$-CROSS constellation is a (rotated) version
of a $36$-AM square constellation with the corner points deleted. Rotation
does not change the energy and so we can get the answer for $32$-CROSS by
finding the total energy for a $36$-QAM square constellation, subtracting
off the energy of the $4$ corner points, and dividing what is left by $32$.


I'm assuming you're using Matlab. The Matlab function and its parameters are:
spectrogram(x,window,noverlap,nfft)

The color indicates the power at a given frequency. Typically, this is given in dB. In your plot, the power is concentrated in the low frequencies a the bottom of the plot. The higher frequencies at the top of the plot have less power.

The value of 10 you used is the width (number of samples) of your windowing function. It takes a vector of n samples and applies weights to the samples such that the ends of the vector taper to zero. This is done to limit leakage (a loss of frequency resolution/introduction of artificial frequencies) during the DFT process. You choose the width of your window according to the frequency resolution you desire (more samples give a finer frequency resolution).

I hope this helps.


The sampling time is the time interval between successive samples, also called the sampling interval or the sampling period, and denoted $T$.

The sampling rate is the number of samples per second. It is the reciprocal of the sampling time, i.e. $1/T$, also called the sampling frequency, and denoted $F_s$.

The frequency axis for the FFT is linked to the number $N$ of points in the DFT and the sampling rate $F_s$. It is defined as $f=k\cdot\frac{F_s}{N}$. With $k$ going up to $N$


ARMA (Auto Regressive Mean Average) Process Representation as AR (Auto Regressive)

Ok,
The answer is as following, any ARMA process which is Minimum Phase system could be represented as MA or AR system.

The MA part is easy by Wold's Decomposition.
Since the filter in Wold's Decomposition is Invertible (Minimum Phase by itself), Shall be called $ H \left( z \right) $, It is easy to build the optimal linear predictor from it, $ P \left( z \right) = 1 - \frac{1}{H \left( z \right)} $, which a linear combination of the past of the process -> AR Model.

It is all coming from the property of the Minimum Phase system being analytic in the unit circle.


Looks like aliasing.

Your 360 Hz is probably not a real component at 360 Hz but it's the 9th harmonic (540 Hz) that aliases back to 360 Hz at a 900Hz sample rate.

That's easy to test: change the sample rate to, 940 Hz and see of the 360Hz moves to 400Hz or stays put. It moves to 400Hz, it's clearly aliasing.


The documentation clearly says (emphasis mine):

  If –s option is specified, the score calculated by
  dynamic time warping, that is, the distance between the test data and the reference data is output and sent to Scorefile.

That means that your decimal number, $0.678382$, is the actual DTW distance between your 2 speech samples.

Now DTW is a distance measurement, and as any distance it only gives you relative information. For example, you could set a treshold and says that under this treshold, the 2 signals correspond to the same word.  If you apply DTW on a pitch contour, you could determine whether the 2 melodies are similar (that's used in query by humming).

Of course, the physical meaning of the distance depends completely on the one of the time series compared : if you compare voltage values, or frequency values, etc.
In your case, you compare normalized volumes mapped to unsigned 8-bit data.

I advise you to learn more about DTW and to try to analyse what the distance actually tells you, and how that relates to your application


How to analyse anti-causal discrete transfer function using matlab?

The transfer function
$$H(z)=-z^{-2} -2z^{-1} +2z +z^2$$
can be written as
$$\begin{align}
H(z)&=z^2z^{-2}\left(-z^{-2} -2z^{-1} +2z +z^2\right)\\
&=z^2\left(-z^{-4} -2z^{-3} +2z^{-1} +1\right)\\
&=H_1(z)H_2(z)
\end{align}$$
 It can now be seen more easily that

$H_1(z)$ has two zeros at $z=0$ (and two poles at $\pm\infty$) .
$H_2(z)$ has four poles at $z=0$ and four zeros: one at $z=1$ and three at $z=-1$. What can be said about pole-zero map of $H(z)$ overall?
Since $|H_1(e^{j\omega})|=1$, we cnclude that $|H(e^{j\omega})|=|H_2(e^{j\omega})|$
Since $\angle H_1(e^{j\omega})=2\omega$, we cnclude that $\angle H(e^{j\omega})=\angle H_2(e^{j\omega})+2\omega$


You want to remove the heart beat signal and keep the "noise". We can solve this problem by using a denoising algorithm, and subtracting the denoised signal from the original signal.

Setting frequency cutoffs for a frequency domain filter can get tricky and turn into a game of whack-a-mole because there's "high frequency" components in the heartbeat blip (due to the sudden rise and fall) and also in the wiggly stuff between a heart beat's T wave and the next heart beat's P wave.

Loosely speaking, the requirements in this de-noising problem are as follows:

Remove the little wiggles
Maintain larger jumps that appear in a heart beat PQRST waveform

This sounds like a great place to apply $\ell_1$ denoising or total variation denoising. The idea is to approximate the given signal $y$ with a signal $x$ such that the derivative of $x$ is "sparse" i.e. it doesn't change too frequently, but when it changes, the change is large. The denoised estimate $\hat x$ is given by:
$$
\hat x = \min_x ||y-x||_2^2 +\lambda\sum|x_i-x_{i-1}|
$$

I used proxTV toolbox in Python to solve this optimization problem.

import scipy.io as sio
import numpy as np
import matplotlib.pyplot as plt
import prox_tv as ptv

mat_struct = sio.loadmat('Signal1.mat')
noisy_signal = mat_struct['x'].T[0]

filtered_signal = ptv.tv1_1d(noisy_signal, 50)

time_vec = np.linspace(0, len(noisy_signal)/1500., len(noisy_signal))

plt.close('all')

fig, ax = plt.subplots(3,1,sharex=True)

ax[0].plot(time_vec,noisy_signal)
ax[0].set_title('noisy signal')

ax[1].plot(time_vec,filtered_signal)
ax[1].set_title('filtered signal')

ax[2].plot(time_vec,noisy_signal - filtered_signal)
ax[2].set_title('noise')
ax[2].set_xlabel('time (s)')

plt.tight_layout()
plt.show(block=False)

And here's the resulting plot:

Of course, there's a different kind of whack-a-mole you'll have to play with this technique: the $\lambda$ parameter which I set to $50$ in the code.


How to derive the stationary Kalman filter predictor?

Your derivations are correct.

$\bar P = P(t|t-1)$ and $K(t) = A \bar K$

Is this your confusion:

Why didn't they have the term $t|t-1$ in the Kalman Gain and Covariance Matrix Expressions?
How can this be "stationary" when your derivation shows that it is time varying?

Bad choice of notation on the book's part

Let's look at the expression:
$ \bar P = A\bar PA^T - A\bar P C^T(C\bar P C^T + R)^{-1}C\bar PA^T + Q $.
The fact that $\bar P $ is a function of itself shows a recursive relationship. In other words, it uses its past values. So, it is NOT the same for all time instants -  It changes at every iteration.

Misunderstanding of the word "stationary".

When the author of the book said "stationary" he/she did not mean that $P$ and $K$ have the same value at all time instants. Instead, the author wanted to emphasise that the expressions for those to values are the same for all statistical realisations.  Stationarity is a statistical concept which means that the Statistics of the system is same all the time. Look at the expressions for $\bar P$ and $\bar K$ again. They depend only on
\

Previous values of themselves
Transition matrices $A$ and $C$ which are deterministic and in your case time-invariant ($A$  and $C$ are the same at all times)
$Q$ and $R$ which are the noise covariance matrices. These 2 matrices describe the statistics of the noises and are the same in all realisations and time instances.

The Kalman gain, $K$, and the state covariance Matrix $P$ will have the same value for all realisations of this random process. (Side Note: None of these 2 terms depend on the measurements, $y$. So they can be computed beforehand.)

Conclusion:

The "time-variant" equations you derived were equivalent to the ones in the book. Besides, the notational differences, there was a slight misunderstanding on your part regarding what changes and what doesn't.


How to set the frequency cut-off for a low pass RC filter?

You probably know that the frequency response of this first order low-pass filter is

$$H(\omega)=\frac{1}{1+j\omega RC}$$

Now compute its magnitude

$$|H(\omega)|=\frac{1}{\sqrt{1+\omega^2R^2C^2}}$$

In order to find the 3-dB cut-off point $\omega_0$, you have to solve $|H(\omega_0)|=1/\sqrt{2}$:

$$\frac{1}{\sqrt{1+\omega_0^2R^2C^2}}=\frac{1}{\sqrt{2}}\Rightarrow
\omega_0=\frac{1}{RC}$$

And since the frequency in Hertz is $f=\omega/2\pi$, you get

$$f_0=\frac{1}{2\pi RC}$$


find the maximum of the DFT of sampled signal

Note that the frequencies of the FFT grid are given by

$$f_k=\frac{kf_s}{N},\qquad k=0,1,.\ldots,N-1\tag{1}$$

where $N$ is the FFT length.

Now observe that the two frequencies of the given signal fall exactly on the grid. So you only have to determine these two frequencies, and then use $(1)$ to figure out the corresponding indices.

EDIT:

It's always good to check your results using some simulation software. The following Octave/Matlab script confirms what you came up with.

phi = 0;    % exact value irrelevant
n = 0:99;
T = 1/25;
y = cos( 10*pi*n*T+phi ) + cos( 20*pi*n*T );
Y = fft(y);
subplot(2,1,1), plot(n,abs(Y),'x'), grid on


Sampling in the time domain is equivalent to periodic repetition in the frequency domain with $1/T_s$ as the period. So formally you'd have something like
$$R_{sampled}(f) = \sum_{n=-\infty}^{\infty}R_z(f-2\cdot n)$$

If you evaluate this, you actually get $R(f) = 2$. It's easier to see if you draw it: take the initial spectrum, shift it by -6,-4,-2,0,2,4,6... and sum it all graphically


How to find orthonormal basis function in the following digital communication problem?

Let's call your $4$ signals $f_0$, $f_1$, $f_2$ & $f_3$.
Just by eyeballing you could have a basis of two function where

b0 = f0, for 0 < t < Ts/2, 0 otherwise and
b1 = f0, for Ts/2 < t < Ts, 0 otherwise and

So you basically have a one first half ramp up in the first half and one second half ramo down. The signals can be expressed simply as

f0 =  b0 + b1
f1 = -b0 - b1
f2 =  b0 - b1
f3 = -b0 + b1

So it's clearly a base. It's also orthogonal since $b_0  b_1$ is zero for all $t$. It's probably not normal but could be easily made so by proper scaling.


How to compute camera pose from Homography matrix?

Even if my answer comes too late for you, maybe other people find this useful. I have the codes for an openCV Pose from Homography. I found the method at this really useful website, euclideanspace.

void cameraPoseFromHomography(const Mat& H, Mat& pose)
{
    pose = Mat::eye(3, 4, CV_64FC1); //3x4 matrix
    float norm1 = (float)norm(H.col(0));
    float norm2 = (float)norm(H.col(1));
    float tnorm = (norm1 + norm2) / 2.0f;

    Mat v1 = H.col(0);
    Mat v2 = pose.col(0);

    cv::normalize(v1, v2); // Normalize the rotation

    v1 = H.col(1);
    v2 = pose.col(1);

    cv::normalize(v1, v2);

    v1 = pose.col(0);
    v2 = pose.col(1);

    Mat v3 = v1.cross(v2);  //Computes the cross-product of v1 and v2
    Mat c2 = pose.col(2);
    v3.copyTo(c2);

    pose.col(3) = H.col(2) / tnorm; //vector t [R|t]
}

////


Convolution of Input signal and Impulse response

In the second case, which caused you the trouble, you should first use the sifting property of the impulse to find out the input signal $x[n]$ and then proceed with convolution.

sifting property is:
$$ f[n] \delta[n-d] = f[d]\delta[n-d]$$

Apply this to simplify your input signal $x[n]$ as follows:

\begin{align}
x[n] &= (\frac{1}{4})^{-n} \left(\delta[n+1] - \delta[n-4] \right) \\
& = (\frac{1}{4})^{-n} \delta[n+1] - (\frac{1}{4})^{-n} \delta[n-4] \\
& = (\frac{1}{4})^{-(-1)} \delta[n+1] - (\frac{1}{4})^{-(4)} \delta[n-4] \\
x[n] & =  \frac{1}{4} \delta[n+1] - 256 ~ \delta[n-4] \\
\end{align}

So, your input sequence is $x[n] = \frac{1}{4} \delta[n+1] - 256 ~\delta[n-4] $

Convolving this with the impulse response of
$$h[n] = \delta[n] - \delta[n-5]$$ is easy, this time by utilizing the shifting property of the impulse which is:

$$ f[n] \star \delta[n-d] = f[n-d] $$

Where $\star$ is th convolution. Then apply distribution of convolution over addition to get your result as:

\begin{align}
y[n] &= x[n] \star h[n] \\
& = x[n] \star (\delta[n] - \delta[n-5] ) \\
&= x[n] - x[n-5] \\
&= (\frac{1}{4} \delta[n+1] - 256 ~\delta[n-4]) - (\frac{1}{4} \delta[n-4] - 256~ \delta[n-9]) \\
&= \frac{1}{4} \delta[n+1] - 256~ \delta[n-4] - \frac{1}{4} \delta[n-4] + 256~ \delta[n-9] \\
y[n] &= \frac{1}{4} \delta[n+1] - (256~ \delta[n-4] + \frac{1}{4}) \delta[n-4] + 256~ \delta[n-9] \\
\end{align}

I hope this helps.


What is the bandwidth of product of two signals?

The bandwidths simply add.

You can break down both signal into their sinusoidal components and pairwise multiply them. For every pair of sines you get the sum and difference frequencies. This is a simple consequence of the multiplication theorems for sine waves. The highest frequency of the product will be the sum of the highest frequencies in the individual signals. This is equivalent to a convolution in the frequency domain.


How to compare many bandpass filters (for their signal propogation delay)?

For a linear phase FIR filter, the group delay will be about the same delay as what you have defined to be a rectengular envelope delay.

For a nonlinear phase IIR or FIR filter (like the Butterworth), the envelope delay is associated with the dominant time constant of its impulse respone $h[n]$.

Indeed what's happening in your plot is that the filter reaches its steady-state response, after the transients decay. And that would happen after a few time constants. For the relation between group delay of the filter and the time constant, it can be said roughly that, short transition bandwidths in the frequence response (sharp filter or narrow band filter with large time constants) yield more aggressive phase changes which also yields larger group delays at the passband frequencies.

A minimum-phase filter would yield the minimum delay (among the same frequency response magnitudes). A minimum-phase filter's impulse response has the fastest decay and begins more aggresively than the other filters; i.e., it's energy is concentrated at the very beginning of its impulse response samples.


Brownian noise is integration of gaussian or uniform white noise?

White noise is noise that has equal (uniform) amplitude across all frequencies. When we say "white" we're talking about the power spectral density (PSD) of the noise.

Saying something like "Gaussian noise" means the statistical properties of any one sample of the noise is distributed Gaussian. You can actually have Cauchy, Poisson, Gaussian etc. distributions that define any one sample of the noise. For any distribution, however, if the power spectral density of the noise is uniform, that noise is white.

So in the integral defining the Weiner process:
$$W(t) = \int_0^t\frac{dW(\tau)}{d\tau}{d\tau}$$

$W(t)$ has samples with Cauchy, Poisson, Gaussian, etc. distributions as long as the power spectral density is uniform, making it white noise.


extract the tree structure from a noisy data set

I guess a very good resource about that is here (Implementation might be tough, but even reading it will give you a lot of insight about the extraction of 3D centerlines)

Also, please take a look at 1 and 2.


'xcorr' function in matlab not working for zero signal

the xcorr function is defined as
$$R_{xy}(m)=
\left\{
 \begin{array}{ll}
  \sum\limits_{n=0}^{N-m-1}x_{n+m}y^*_n  & \mbox{if } m \geq 0 \\
  R^*_{yx}(-m)& \mbox{if } m < 0
 \end{array}
\right.$$
The output is given as $$c(m)=R_{xy}(m-N)$$ where N is the length of the signal by the above definition the output for zero signal is zero


Inverse of Wavelet Transforms - Background and Noise Removal

I ended up finding a package which provides an inverse wavelet transform function called mlpy.  The function is mlpy.wavelet.uwt.  This is the snippet which may interest people if they are trying to do noise or background removal:

# Make 2**n amount of data
new_y, bool_y = wave.pad(xy['y'])
orig_mask = np.where(bool_y==True)

# wavelet transform parameters
levels = 8
wf = 'h'
k = 2

# Remove Noise first
# Wave transform
wt = wave.uwt(new_y, wf, k, levels)
# Matrix of the difference between each wavelet level and the original data
diff_array = np.array([(wave.iuwt(wt[i:i+1], wf, k)-new_y) for i in range(len(wt))])
# Index of the level which is most similar to original data (to obtain smoothed data)
indx = np.argmin(np.sum(diff_array**2, axis=1))
# Use the wavelet levels around this region
noise_wt = wt[indx:indx+1]
# smoothed data in 2^n length
new_y = wave.iuwt(noise_wt, wf, k)

# Background Removal
error = 10000
errdiff = 100
i = -1
iter_y_dict = {0:np.copy(new_y)}
bkg_approx_dict = {0:np.array([])}
while abs(errdiff)>=1*10**-24:
    i += 1
    # Wave transform
    wt = wave.uwt(iter_y_dict[i], wf, k, levels)

    # Assume last slice is lowest frequency (background approximation)
    bkg_wt = wt[-3:-1]
    bkg_approx_dict[i] = wave.iuwt(bkg_wt, wf, k)

    # Get the error
    errdiff = error - sum(iter_y_dict[i] - bkg_approx_dict[i])**2
    error = sum(iter_y_dict[i] - bkg_approx_dict[i])**2

    # Make every peak higher than bkg_wt
    diff = (new_y - bkg_approx_dict[i])
    peak_idxs_to_remove = np.where(diff>0.)[0]
    iter_y_dict[i+1] = np.copy(new_y)
    iter_y_dict[i+1][peak_idxs_to_remove] = np.copy(bkg_approx_dict[i])[peak_idxs_to_remove]

# new data without noise and background
new_y = new_y[orig_mask]
bkg_approx = bkg_approx_dict[len(bkg_approx_dict.keys())-1][orig_mask]
new_data = diff[orig_mask]

Here a complete compilable script which can produce output:

import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
import mlpy.wavelet as wave

# Make some random data with peaks and noise
############################################################
def gen_data():
    def make_peaks(x):
        bkg_peaks = np.array(np.zeros(len(x)))
        desired_peaks = np.array(np.zeros(len(x)))
        # Make peaks which contain the data desired
        # (Mid range/frequency peaks)
        for i in range(0,10):
            center = x[-1] * np.random.random() - x[0]
            amp = 100 * np.random.random() + 10
            width = 10 * np.random.random() + 5
            desired_peaks += amp * np.e**(-(x-center)**2/(2*width**2))
        # Also make background peaks (not desired)
        for i in range(0,3):
            center = x[-1] * np.random.random() - x[0]
            amp = 80 * np.random.random() + 10
            width = 100 * np.random.random() + 100
            bkg_peaks += amp * np.e**(-(x-center)**2/(2*width**2))
        return bkg_peaks, desired_peaks

    # make x axis
    x = np.array(range(0, 1000))
    bkg_peaks, desired_peaks = make_peaks(x)
    avg_noise_level = 30
    std_dev_noise = 10
    size = len(x)
    scattering_noise_amp = 100
    scat_center = 100
    scat_width = 15
    scat_std_dev_noise = 100
    y_scattering_noise = np.random.normal(scattering_noise_amp, scat_std_dev_noise, size) * np.e**(-(x-scat_center)**2/(2*scat_width**2))
    y_noise = np.random.normal(avg_noise_level, std_dev_noise, size) + y_scattering_noise
    y = bkg_peaks + desired_peaks + y_noise
    xy = np.array( zip(x,y), dtype=[('x',float), ('y',float)])
    return xy
# Random data Generated
#############################################################

xy = gen_data()

# Make 2**n amount of data
new_y, bool_y = wave.pad(xy['y'])
orig_mask = np.where(bool_y==True)

# wavelet transform parameters
levels = 8
wf = 'h'
k = 2

# Remove Noise first
# Wave transform
wt = wave.uwt(new_y, wf, k, levels)
# Matrix of the difference between each wavelet level and the original data
diff_array = np.array([(wave.iuwt(wt[i:i+1], wf, k)-new_y) for i in range(len(wt))])
# Index of the level which is most similar to original data (to obtain smoothed data)
indx = np.argmin(np.sum(diff_array**2, axis=1))
# Use the wavelet levels around this region
noise_wt = wt[indx:indx+1]
# smoothed data in 2^n length
new_y = wave.iuwt(noise_wt, wf, k)

# Background Removal
error = 10000
errdiff = 100
i = -1
iter_y_dict = {0:np.copy(new_y)}
bkg_approx_dict = {0:np.array([])}
while abs(errdiff)>=1*10**-24:
    i += 1
    # Wave transform
    wt = wave.uwt(iter_y_dict[i], wf, k, levels)

    # Assume last slice is lowest frequency (background approximation)
    bkg_wt = wt[-3:-1]
    bkg_approx_dict[i] = wave.iuwt(bkg_wt, wf, k)

    # Get the error
    errdiff = error - sum(iter_y_dict[i] - bkg_approx_dict[i])**2
    error = sum(iter_y_dict[i] - bkg_approx_dict[i])**2

    # Make every peak higher than bkg_wt
    diff = (new_y - bkg_approx_dict[i])
    peak_idxs_to_remove = np.where(diff>0.)[0]
    iter_y_dict[i+1] = np.copy(new_y)
    iter_y_dict[i+1][peak_idxs_to_remove] = np.copy(bkg_approx_dict[i])[peak_idxs_to_remove]

# new data without noise and background
new_y = new_y[orig_mask]
bkg_approx = bkg_approx_dict[len(bkg_approx_dict.keys())-1][orig_mask]
new_data = diff[orig_mask]

##############################################################
# plot the data and results
fig = plt.figure()

ax_raw_data = fig.add_subplot(121)
ax_WT = fig.add_subplot(122)

ax_raw_data.plot(xy['x'], xy['y'], 'g')
for bkg in bkg_approx_dict.values():
    ax_raw_data.plot(xy['x'], bkg[orig_mask], 'k')

ax_WT.plot(xy['x'], new_data, 'y')

fig.tight_layout()
plt.show()

And here is the output I am getting now:

As you can see, there is still a problem with the background removal (it shifts to the right after each iteration), but it is a different question which I will address here.


Help me understand the stages involved in filtering a signal using Discrete Fourier Transform

A purely real time domain signal has a frequency transform such that positive and negative frequencies are complex conjugated of one another.  If you are scaling the positive bin, but not the negative, you are implementing a complex filter kernel, resulting in a complex output.

Side note, you can’t just change the magnitude of a frequency bin in a frame and do an inverse transform in real time processing.  The phases of adjacent bins won’t line up.  This is a more advanced topic beyond the scope of your question.


Are capacity and spectral efficiency identical for a dicrete-time digital signal?

The capacity is a property of a channel. A channel can transport no more than a certain amount of information per channel use. When we use a fixed modulation and coding scheme ("discrete-time signal"), this scheme will be able to support a certain amount of information per channel use (which for the channel to support it must be $\leq$ its capacity). It is thus strictly speaking incorrect to speak of the capacity of a certain modulation and coding scheme.
For a given modulation and coding scheme such as a QAM, you can of course compute its achievable data rate in bits per channel use (i.e., per unit bandwidth). You can also interpret this as a spectral efficiency, since it measures how efficient the scheme uses bandwidth (how many Mbit/s do I get for each MHz of bandwidth I invest?). And the results you get may be reminiscent of the capacity since under certain conditions, the schemes can get somewhat close to it.
It is still not the same. Capacity is a property of a channel, spectral efficiency is a property of a transmission scheme (modulation and coding scheme).


state space formulation of a sinusoidal system

I make you an example about how to write state-space equations using a dynamic second order equation to let you understand the procedure.

The starting equation is: $$M a + C v + K  s = f $$
$C$ is damping, $K$ is stiffness and $M$ is the mass.

Then build your state vector $x =\pmatrix{s\\ v}$

Now write the equation like $ a = 1/M (f - C v - K s)$

$\dot{x} = \pmatrix{v\\ a}$
is the derivative of the state vector.

In the state space:

$\dot{x} = \pmatrix{0& 1, \\-K/M& -C/M} x + \pmatrix{0 \\ 1/M} f$

If you want the displacement as output

$y = \pmatrix{1\\ 0} x $


An obvious choice would be $\frac{d_1}{d_1+d_2+\epsilon }$, where $\epsilon$ is a small value to prevent division by zero.

Using your images, I get:

top = ColorNegate[Binarize[ImageApply[StandardDeviation,
           Import["http://i.stack.imgur.com/aFNfB.jpg"]]]];
bottom = ColorNegate[Binarize[ImageApply[StandardDeviation,
           Import["http://i.stack.imgur.com/Kipvo.jpg"]]]];
tDist = ImageData[DistanceTransform[top]];
bDist = ImageData[DistanceTransform[bottom]];

\[Epsilon] = 10^-10;
Image[tDist/(tDist + bDist + \[Epsilon])] // ImageAdjust


Why does this not work: Alt. way to compute the period of a discrete periodic signal that is a sum of 2 complex exponentials

The flaw in your reasoning is that you only consider the cosine term. The cosine term has a period of $48$. But there's also the complex exponential, which is not just a magnitude but a function of $n$. You can show that that exponential term is periodic with a period of $48$. However, you cannot conclude that the multiplication of the two terms must also have a period of $48$.

Multiplying the two terms results in a period of $24$ because

$$e^{j\frac{8.5\pi}{12}n}=-e^{j\frac{8.5\pi}{12}(n+24)}$$

and

$$\cos\left(\frac{0.5\pi n}{12}\right)=-\cos\left(\frac{0.5\pi (n+24)}{12}\right)$$

so the sign changes cancel in the product.


Estimate the Discrete Fourier Transform / Series of a Signal with Missing Samples

Given $ \left\{ x \left[ n \right] \right\}_{n \in M} $ where $ M $ is the set of indices given for the samples of $ x \left[ n \right] $.

The trivial solution (Which it would be great to have a faster more efficient solution is what I'm looking for) would be:

$$ \arg \min_{y} \frac{1}{2} \left\| \hat{F}^{T} y - x \right\|_{2}^{2} $$

Where $ \hat{F} $ is formed by subset of columns of the DFT Matrix $ F $ matching the given indices of the samples, $ x $ is the vector of the given samples and $ y $ is the vector of the estimated DFT of the full data of $ x \left[ n \right] $.

The solution is then given by the Pseudo Inverse (Least Squares Solution):

$$ y = { ( \hat{F} \hat{F}^{T} ) }^{-1} \hat{F} x $$

In practice, the matrix will be very poorly conditioned hence solution must be generated using the LS Solution using the SVD.

A sample code is shared on GitHub Repository.

Result of the code:


What is the purpose of the recurrence relation in low and high pass (audio) filters?

An IIR filter is a recursive one. It means, current value of output $y[n]$ is computed based on its previous value(s).
Consider a simple 1st order IIR filter definition:
$$y[n] = a \cdot y[n-1]+ b \cdot x[n]$$
$a$ and $b$ are the filter coefficients.
Now assume you have a long block of N samples input data. If you want to process that block of input data with this filter, you can do it like in this matlab code :
N = 1024;       % number of input samples
x = randn(N,1); % some random input.
a = 0.5;        % IIR filter coefficients
b = 0.9;

yi = 0;         % initial value for output
y = zeros(N,1); % matrix (array) to hold processed output samples

y(1) = a*yi + b*x(1)   % assign first output sample based on y(-1)
for n=2:N              % assign rest
    y(n) = a*y(n-1) + b*x(n);
end


Rake Receiver vs OFDM and multi-path channel

You can basically apply the rake receiver whenever.

However, OFDM was invented especially to make that unnecessary.

The point of a rake receiver is to demodulate the different paths independently and sum up the observed information. It is a time-domain decomposition of the channel impulse response, if you will.

OFDM considers the effects of the multipath channel, but in frequency domain, and converts the complicated time-domain equalization problem (which the Rake, partially, solves) to a frequency-domain problem, where the channel properties are independent between OFDM subcarriers and trivial to correct. It's the frequency-domain decomposition of the channel, if you will.

So, if you do OFDM, and your channel delay spread is within what your OFDM system can deal with, adding a rake in before seems redundant.

Your "pick the stronges rake path" is simply like a bad OFDM synchronizer. So, this is worse than just doing the very standard Schmidl&Cox synchronization for packetized OFDM, or to fixed-lag correlation methods used in streaming OFDM applications.

  regarding how can I get the BER, I can do it similarly to this paper journals.sagepub.com/doi/pdf/10.1155/2016/3137014 .. But anyway, don't care of that method, let's consider the method of selecting the strongest signals. It means, "the signal can reach the receiver via multiple distinct pathways. The rake receiver is used to correct this effect, selecting the stronger signals.Each of these multipath components in rake receiver is called finger" so then, we can demodulate the different taps in different demodulators

Yeah, no. An OFDM synchronization would already do that for you. Rake and OFDM solve the same problem: Multipath == frequency selective channels. But OFDM inherently solves the problems of single carrier systems; applying Rake before doing OFDM is IMHO a dead end.

You need to read up on both when to use rake, and when to use OFDM. Recent papers aren't the way here; they are either low-quality (explaining concepts that have been standard for decades), or extremely special in application. Get a good communications textbook!


What's Logic Behind the Construction of Sobel's Filter in Image Processing?

A first rationale is to be very short, as there was a time when computing on images was expensive. Then, a contour or an edge often present a fast variation in image intensities, that can be enhanced by derivatives. Sobel filters emulate such derivatives in one direction, and slightly average pixels in the complementary direction, to smooth small variations or noise.

One direction implements the shortest possible centered 1D discrete derivative:

$$\begin{bmatrix}
   -1 &0 &1
\end{bmatrix}
$$
to detect variations across lines, the other the shortest non-trivial Pascal/Gaussian smoothing
$$\begin{bmatrix}
   1&2&1
\end{bmatrix}^T
$$
to smooth along columns, resulting in, for instance:
$$
\begin{bmatrix}
   1&2&1
\end{bmatrix}^T\cdot
\begin{bmatrix}
   -1 &0 &1
\end{bmatrix}
$$
or
$$
\begin{bmatrix}
   -1 &0 &1 \\
   -2 &0 &2 \\
   -1 &0 &1 \\
\end{bmatrix}
$$

As you can see, this only involves dyadic numbers, so it can be implemented with adds and binary shifts.

Of course, the 3-point derivative often has an additional $1/2$ factor: $$\begin{bmatrix}
   -1/2 &0 &1/2
\end{bmatrix}
$$
to get the appropriate scale factor, and the Pascal smoother has a $1/4$ factor to have its coefficients sum to one
$$\begin{bmatrix}
   1/4&1/2&1/4
\end{bmatrix}
$$
but the resulting global scaling of $1/2\times 1/4$ does not change the edge detection power for such linear filters.


More Simultaneous Dynamic Range with Fixed Point or Floating Point?

For the same total bit width, a fixed-point representation has a higher instantaneous dynamic range than a floating-point representation.
 For a value to be simultaneously including information about a strong and a weak signal (where the information is being carried by the sum of the strong and weak signal components), the weak signal has to have a value at least as large as the stronger signal's numerical error floor, given by the format's epsilon. In a floating-point representation, every bit allocated to the exponent similarly reduces the number of bits allocated to the significand, and thus increases the format's epsilon. This correspondingly decreases the ratio between the strong and weak signal components that can be included in a single value.
 For example in a standard 32-bit IEEE-754 single-precision floating point representation, 23 bits are allocated to the significand. As a result the smallest number that can be added to a given number $x$ is of the order of $x/2^{23}$ or approximately $8.6\times 10^6$ times smaller. Note that this is consistent with the C programming language's FLT_EPSILON which is approximately $1.19\times 10^{-7}$. In comparison, using a 32-bit fixed-point representation the smallest number that can be added to a large number using near full scale (e.g. 0xffffff00) is of the order of a single unit (e.g. 0x00000001) which is approximately $2^{32} \approx 4.3\times10^9$ times smaller.


Aliasing correspondance with respect to DFT

If you do spectral analysis is typically better to use WAY more points then 10 in the time domain
Look at your y axis and make sure your graphs are properly scaled. $f=10$ aliases indeed back to $f=0$. All you see is numerical noise and your y-scale is $10^{-15}$. This would be different if you choose a a different phase, cosine instead of sine.


Expanding piecewise polynomial using Daubechies wavelet

According to your formula, you also apply soft-thresholding to the approximation coefficients, which is not standard. Aside, your operator $W$ does not seem to specify the number of wavelet levels used. Finally, your class of signals does not seem to address the regularity at piecewise junctions.

I believe in this case very unlikely that in a discrete implementation, without further relation between $m$ and $n$, you can find, theoretically,  a best wavelet in all cases, because $\|\hat{y}\|_0$ is a quite sensitive index (and it is not a norm).  Of course, a Daubechies with $d$ vanishing moments would seem appropriate.

But since the DWT is quite fast, in your context, you could just find a "generally best wavelet",  by simulating random signals, and iterating over some levels, and each Daubechies wavelet with moments in $[d-2,\ldots,d+2]$ for instance.


Is it numerically more stable to implement filtering as multiplication or convolution?

With convolution, you are not going to run into any stability problems, because there is no recursive filtering, so you are not going to accumulate any errors. In other words, the system is all zeros, no poles. I have heard anecdotally, but not verified for myself, that FFT-based convolution has lower error than time-domain convolution, simply because it has O(n log n) arithmetic operations rather than O(n^2).

Typically, as far as I am aware, Butterworth filters are implemented as recursive (IIR) filters, so that is a different topic. IIR filters have poles as well as zeros, so there can be stability issues in practice. Also, for IIR filters, FFT-based methods are not an option, but on the upside, IIR filters tend to be very low order.

As far as stability issues with IIR filters, they tend to have problems at higher orders -- I'll just throw out a number and say roughly 6th order is pushing it. Instead, they are typically implemented as cascaded biquads (2nd order filter sections). For your 5th order filter, write it as a z-domain transfer function (it'll be a 5th degree rational function), and then factor it into its 5 poles and 5 zeros. Collect the complex conjugates, and you'll have two biquads and one first-order filter. In general, stability problems tend to crop up as the poles get closer to the unit circle.

There can also be issues with noise and limit cycles in IIR filters, so there are different filter topologies (i.e. direct form I, direct form II) that have different numeric properties, but I wouldn't overthink this point -- just use double-precision and it'll almost certainly be good enough.


Are there any cases where the negative frequencies of an FFT are used?

A strictly real vector in one domain implies a complex conjugate mirrored vector in the other domain, and vice versa, for any FFT.

The "negative" frequency results of an FFT are very commonly used in SDR (Software Defined Radio), where the signal of interest might consist of only negative frequencies (LSB or lower sideband SSB is one obvious case) in a baseband IQ stream (after quadrature downconversion from strictly real higher frequency RF).  Or the SDR software may display a wideband spectrum showing signals both above and below the IQ baseband center frequency, which are displayed from an FFTs negative frequency bins thru positive frequency bins.

Some USB SDR devices (Funcube, et.al.) communicate the complex IQ data stream as stereo audio inputs to the computer, so one might actually be able to listen to the complex IQ data as a stereo left and right pair of audio (after appropriate sample rate conversion, if needed).


Numericaly computing an integral in MATLAB

Use integral function.

q = integral(g,xmin,xmax)

For example $g(t) = e^{-t}$

$$q = \int_0^{\infty}e^{-t}\mathrm{d}t = 1$$

g = @(t) exp(-t);
q = integral(g,0,Inf)

More information
Matlab numerical integration


Estimate frequency and peak value of a signal's fundamental

There are a few Matlab implementations of various frequency estimators here.  However, they will probably be "tricked" by the higher harmonics.

That you know a range, perhaps the best technique is one of Eric Jacobsen's estimators. Have a look at his page.

In essence, Eric's estimators use the FFT bin at $X(k)$, and on either side of ($X(k+1), X(k-1)$) the peak.  If you choose your FFT length to match your known range (so that $k \rightarrow 3.45$ and $k+1 \rightarrow 3.50$ and $k-1 \rightarrow 3.4$, then you should be good to go.


Removing noisy lines from image - opencv - python

first of all remember that there is no single solution for all kind of noise and all kind of images. that being said i can think of two solution.
first is using Otsu thresholding:

ret,thresh_img = cv2.threshold(img, 0, 255, cv2.THRESH_BINARY_INV|cv2.THRESH_OTSU)

this will try to guess a good threshold for the image being used.

the other solution would be using "close" morphology transform after thresholding. it will first dilate and then erode the image and using a good kernel( it is also called structuring element ) you can remove the line(although you may remove some of the useful pixels too!!)

kernel = cv2.getStructuringElement(cv2.MORPH_RECT, (4,8))
morph_img = cv2.morphologyEx(thresh_img, cv2.MORPH_CLOSE, kernel)

above is the result of using both otsu threshold and morphology close


Reversing the Order of Operators for Edge Detection?

In the classic framework both the Smoothing and the Difference Filter are applied using Convolution.
Since it is done using convolution it implies the operation is Linear Spatially Invariant (LSI).
LSI operators can be applied in any order and the result will be the same. This is also a result of the commutativity property of the convolution operator.
Let's say $ f $ is the smoothing operator and $ g $ is the difference operator then:
$$ y = f \ast g \ast x = g \ast f \ast x $$
In practice, a more computational efficient way is define $ h = f \ast g = g \ast f $.
See Difference of Gaussians.


When talking about image deconvolution, it is referring to the process of correcting an image for distortion caused by the optical properties of the camera or imaging system. This distortion is characterized by the point spread function of the optical system. The point spread function describes the response of an imaging system to a point source. For a non point source ,the image is a superposition of the point spread function of multiple point sources representing the non point source.

There are many techniques of image deconvolution and references, too many to list in this format. But for examples of applications using linear algebra, you may want to begin by looking for references on image deconvolution using
Wiener filters, or the Richardson-Lucy method. Both require knowledge of the point spread function of  the system , but both involve linear algebra.


What is the duration of the spreading sequence

When you multiply the real and imaginary parts by $\mathbf{W}_1$ and $\mathbf{W}_2$, respectively, the resulting signal is no longer $S$. It will be a matrix.

Regarding your questions:

1- Yes, the signal after spreading has the same duration.

2- One of the disadvantages of Walsh-Hadamard codes is that they don't have clear auto-correlation peak, which makes synchronization at the receiver side challenging. Also, they must have a dimension of $2^k$, i.e., you cannot have arbitrary dimension to support an arbitrary number of users. In your example, to make the system efficient, you need 4 users, otherwise, one user will occupy 4 times the minimum required bandwidth, which isn't spectrally efficient.

3- Yes. This is called direct sequence (DS) CDMA.


2 points and instantaneous phase sufficient to estimate amplitude and frequency of sinusoid?

An additional requirement is that the two points not be aliased in terms of their equivalent sampling rate relative to the frequency of the sinusoid.  The distance between the two points can't be a equal or greater than half the period of the sinusoid, or the solution could be ambiguous.


How can i compute perimeter and area of a skeletonized circle?

It is not clear if the skeletonized image can be reversed to a contour image. Skeletonization is often done by calculating the distance transform and then thresholding the result or performing a non-maximum suppression. At this point you have already lost too much information and the features, you are asking for cannot be calculated.

If you have the distance transform, you can approximate those values:

You get the radius ($r$) of the circle from the height of peak which represents the circle (the skeleton of an ideal circle would be a dot)
Perimeter = $2\pi r$
Area = $\pi r^2$
Convex hull $\rightarrow$ draw a circle with radius $r$ at the position where the peak is


Defining the observation matrix for a Linear Kalman Filter

The first question you link to states that it only measures acceleration and the answer suggests to add the acceleration to the state vector as well, $\begin{bmatrix}x^\top&\dot{x}^\top&\ddot{x}^\top\end{bmatrix}^\top$. This is why the observation matrix becomes $\begin{bmatrix}0&0&I\end{bmatrix}$.

So your reasoning for observation matrix is correct. I only have one remark, namely it is given that the drone moves in 3D space, so the observation matrix should actually be

$$
\begin{bmatrix}
I&0\\0&I
\end{bmatrix}
$$

with $I$ a three by three identity matrix. So effectively the observation matrix would be a six by six identity matrix.


Using Local Binary Pattern for finding Homography between two Image

Yes. LBP can be used to find homography between two images. LBPs extract the local difference between nearby pixels with respect to center pixel and represent in 8-bit. You can analyse this 8-bit and extract just the corners which can be used as feature points. If you have an option to use other algorithm, I suggest you to go with SIFT .


Given the 4 corners of two QR codes on a photographed sheet of paper, can we get a deskewed rectangular image?

The solution is to use an "homography" / aka "projective transformation" (see this PDF, page 16).

Here is a working code showing how to do it with Python + OpenCV.


how to eliminate interference in the same center frequency

First of all, the CC1101 supports a number of other frequencies, so as a trivial solution, you can switch to using those.

The other thing that you can do at the RF module level is to "shift" your data out of the problematic region. Now, I am assuming that the other module is on 433MHz and it is not a similar module with similar operation. Let's assume that the interfering module is different than the one you are using.
What you could do, is shift to 4-FSK with the highest DEVIATION_M (refer to the datasheet) but use it as 2-FSK. What you would be trying to do is use the two symbols that cause the largest deviation to the centre carrier frequency to avoid your data being contaminated by the interference. Usually, these two symbols would be 00 and 11 (notice these pairs). This is inefficient of course, but it might work as a quick workaround for you.

But what might be happening is that the other RF module is also behaving in the same way. It might also be doing Frequency Shift Keying, in which case, none of the frequencies your RF module might try to hop to is "safe" from noise.

In this case, you have to work with high levels of noise.

The number one way to guard against this kind of interference is with the use of coding.

This can take many forms and have various complexities. One of the simplest "encoding" schemes is a Repetition Code coupled with simple Majority Decoding scheme. You might notice above that your effective B-FSK used the symbols 00 and 11 which, effectively, is a bit like using a repetition code.

One of the issues you might face in this case is that you do not have access to the internals of the RF module. In other words, the only thing you would be doing here is repeat bits $N$ times and at the reception end, receive a frame and decide what bit was sent by assigning the frame to the majority of bits received (either $0$ or $1$).

Obviously the question now is, where does the frame start and where does it end?

To deal with this, you are now going to have to embed your data stream on some form of protocol that specifies that your actual message starts with something like 0000000000111111111100001111001101 followed by your encoded bits.

The frame start sequence has equal chances of being contaminated with noise as the rest of the data, but the idea is to be looking in to your incoming data stream for a sequence that looks as close as possible to the frame start sequence.

And thus, progressively, we are going closer and closer to convolutional codes where you would be inserting redundant 0s and 1s in a predictable way in your data stream to be able to decode the most likely sequence to have been sent. That is a bit more involved at the decoding side but convolutional codes are more powerful than plain repetition codes.

Hope this helps.


The opposite of blind-signal separation is not "non-blind signal separation". The thing is, depending on what you know on the signals you will use different algorithms.

For example, if you know that there is a low frequency signal mixed with a higher frequency one, you can use low/high-pass filtering to get the two signals.

There's no real name for what you're looking for because every case is different. The generic way is more or less filtering the signals with the information you know. But the filtering is always different.


Reconstruction of audio signal from Spectrogram: Implementation

Here is an example in Matlab that I found some years ago. It has been turned to Python by one of us, apparently (STFT ISTFT Matlab Python).

function d = stft(x, f, w, h)
% D = stft(X, F, W, H)                            Short-time Fourier transform.
%   Returns some frames of short-term Fourier transform of x.  Each
%   column of the result is one F-point fft; each successive frame is
%   offset by H points until X is exhausted.  Data is hamm-windowed
%   at W pts..
%   See also 'istft.m'.
% dpwe 1994may05.  Uses built-in 'fft'
% $Header: /homes/dpwe/public_html/resources/matlab/RCS/stft.m,v 1.1 2002/02/13 16:15:55 dpwe Exp $

s = length(x);

if rem(w, 2) == 0   % force window to be odd-len
  w = w + 1;
end

halflen = (w-1)/2;
halff = f/2;   % midpoint of win
acthalflen = min(halff, halflen);

halfwin = 0.5 * ( 1 + cos( pi * (0:halflen)/halflen));
win = zeros(1, f);
win((halff+1):(halff+acthalflen)) = halfwin(1:acthalflen);
win((halff+1):-1:(halff-acthalflen+2)) = halfwin(1:acthalflen);

c = 1;

% pre-allocate output array
d = zeros((1+f/2),1+fix((s-f)/h));

for b = 0:h:(s-f)
  u = win.*x((b+1):(b+f));
  t = fft(u);
  d(:,c) = t(1:(1+f/2))';
  c = c+1;
end;

and the inversion:

function x = istft(d, ftsize, w, h)
% X = istft(D, F, W, H)                   Inverse short-time Fourier transform.
%   Performs overlap-add resynthesis from the short-time Fourier transform
%   data in D.  Each column of D is taken as the result of an F-point
%   fft; each successive frame was offset by H points. Data is
%   hamm-windowed at W pts..
% dpwe 1994may24.  Uses built-in 'ifft' etc.
% $Header: /homes/dpwe/public_html/resources/matlab/RCS/istft.m,v 1.1 2002/02/13 16:16:12 dpwe Exp $

s = size(d);
%if s(1) != (ftsize/2)+1
%  error('number of rows should be fftsize/2+1')
%end

cols = s(2);
xlen = ftsize + cols * (h);
x = zeros(1,xlen);

if rem(w, 2) == 0   % force window to be odd-len
  w = w + 1;
end

win = zeros(1, ftsize);

halff = ftsize/2;   % midpoint of win
halflen = (w-1)/2;
acthalflen = min(halff, halflen);

halfwin = 0.5 * ( 1 + cos( pi * (0:halflen)/halflen));
win((halff+1):(halff+acthalflen)) = halfwin(1:acthalflen);
win((halff+1):-1:(halff-acthalflen+2)) = halfwin(1:acthalflen);

for b = 0:h:(h*(cols-1))
  ft = d(:,1+b/h)';
  ft = [ft, conj(ft([((ftsize/2)):-1:2]))];
  px = real(ifft(ft));
  x((b+1):(b+ftsize)) = x((b+1):(b+ftsize))+px.*win;
end;

EDIT: if you want to start from STFT magnitude, then you run into the more complex problem of phase retrieval. A few paper to start from:

Signal reconstruction from STFT magnitude: a state of the art, DAFx 2011
Non-convex phase retrieval from STFT measurements
STFT phase retrieval: Uniqueness guarantees and recovery algorithms


Magnitude Spectrum, different magnitudes same amplitude

To add a bit to the previous answers: the FFT result is "binned", not continuous.  A plot of an FFT result often fools people by being smooth for large enough FFT lengths.  But the result is really more like a bar chart for a sequence of bins, rather than a line graph.

And the bins are not rectangular in shape (in terms of frequency response).  The nearer a spectral frequency peak is to the center of an FFT bin, the more accurate the magnitude result.  You can get a list of your FFT's bin centers by multiplying sample_rate/length_of_FFT by integers from zero up to half the length.  If a frequency peak is not at the exact center of an FFT result bin, the magnitude gets split among multiple bins, which takes height away from the closest peak bin.

In your example, some of your frequencies were closer to bin centers than others.


Amplitude loss during resampling with IIR

Try adding a gain of 16 before your low pass filter, or equivalently using a low pass filter with a passband gain of 16 instead of 1.


Steps for recognizing hand written signs

That's what you would do with images, not so much with handwriting, which needs more preprocessing and robust features because people don't write exactly the same way twice. See for example:

Improving Offline Handwritten Text Recognition with Hybrid HMM/ANN Models
Adaptive Membership Functions for Handwritten Character Recognition by Voronoi-Based Image Zoning
A visual approach to sketched symbol recognition
Feature extraction and classifier combination for image-based sketch recognition
Sketched symbol recognition with auto-completion
HBF49 feature set: A first unified baseline for online symbol recognition
Neural network-based symbol recognition using a few labeled samples


Using this formula for the percentage overshoot, the damping ratio is given by

$$\zeta=\sqrt{\frac{\log(0.4)^2}{\pi^2+\log(0.4)^2}}=0.28$$

which is more or less the same value that you obtained. Consequently, you have an underdamped system. The peak frequency $\omega_{max}$ of the frequency response is related to the damping ratio and to the natural frequency in the following way (link):

$$\omega_{max}=\omega_n\sqrt{1-2\zeta^2}$$

With $\zeta=0.28$ and $f_{max}=500\text{ Hz}$ this results in a natural frequency $f_n=544.51\text{ Hz}$. Apart from a scaling factor, a second order system is completely characterized by $\zeta$ and $\omega_n=2\pi f_n$.


Matched Filters for attenuation detection?

Your interpretation of the matched filter for radar seems almost right, just to be sure, let me reiterate a bit: Think of it this way: we're sending a pulse function $p(t)$. If a target in distance $d$ reflects the signal, we receive $\gamma \cdot p(t-\tau)$, where $\tau = 2d/c$ is the total travelled delay.

The problem is: how to sample such a signal? Let's say your pulse was rectangular, you would need a lot more samples than the pulse width to know its extent and therefore, the proper delay. This is where the matched filter helps.

Applying matched filtering to your received signal, you filter it with $p(-t)$. This gives a signal comprised of the autocorrelation function $\varphi_p(t) = p(t) * p(-t)$, i.e., we obtain $\gamma \cdot \varphi_p(t-\tau)$, more or less. Now, the ACF always has its maximum at $t=0$, which greatly helps for the sampling. Also sampling the ACF at zero gives the optimal SNR: you gather the whole energy of the pulse into the sampling point. It's also called pulse compression for this reason.

Now, back to your question: You can extract both from the received signal, the delay $\tau$ and the reflection coefficient $\gamma$. In fact, this is often done. People call it RCS (radar cross section). Note that it is not so easy to infer material properties from $\gamma$ as $\gamma$ is typically highly direction dependent: the geometry of the object has a lot of effect on how strong the reflected wave is.


Is FIR decimation the fastest possible method for integer downsampling?

Another way to design a FIR filter involves matching the impulse response to an analog prototype's impulse response.  The idea being that  an IIR can often be truncated to FIR. The FIR will not be symmetric  and not have linear phase but  the group delay can be much smaller.
I believe that Friedlander has a paper  on how to do the approximation. There was a Matlab Signal Processing Toolkit function that did this but I haven't had much cause to use it.
Comet also can do nonsymetric FIR designs.

I also suggest that you try using Parks Mclelan if it can give you a shorter filter than using the window method


Autocorrelation - Stochastic vs deterministic processes

You are correct that these two definitions are for different processes.

$R_x\left(\tau\right) = \frac{1}{\sigma^2} \mathbb{E}\left[\left(X_t - \mu\right)\left(X_{t+\tau} - \mu\right)\right]$

is defined for a stochastic process $X$ distributed according to a law $\nu$.

Whilst for a deterministic, finite-energy signal, $x\left(t\right)$ the autocorrelation is:

$R_x\left(\tau\right) = \int_{-\infty}^\infty x\left(\tau+t\right)x\left(\tau\right)dt$

The Wiener-Khinchin theorem will hold for any square-integrable deterministic $x\left(t\right)$ that is:

$\int_{-\infty}^\infty |x\left(\omega\right)|^2 d\omega < \infty$

so that the Fouier integral exists and converges in the $L^2$ sense (i.e. signals with finite energy).

The two definitions are different because the integrals used in the definitions are different. The deterministic case implicitly assumes a Riemann integral - that the process is nice in some way - it is continuous, and if differentiable the derivatives are bounded etc. However for random processes, $dt$ is difficult to define - the process could jump from $-\infty$ to $\infty$ over any inteval I could possibly define, so the Riemann limit isn't appropriate.

Instead we use the Lebesgue integral, where we can define the expectation of a random process: all that's required is that the range is smooth (for example $\mathbb{R}$), and we can still take the expectation with respect to some measure, modulo some restrictions on the random process (so that the integral converges in a monotone or dominated sense).

Integration of Random Process is a tricky business, and the definitions are written differently to keep people mindful of what they are working with.

One other note - the deterministic case isn't so interesting, as we can take the FT of a deterministic signal: the whole point of the WK theorem is for random processes - where the Fourier integral may not exist. I find it a particularly satisfying theorem, and the work to understand it definitely worthwhile.


Using the wav version of the file you linked to I performed an fft in matlab and got:

First 3 terms:

 left channel       right channel
-1.1594 + 0.0000i  -1.8450 + 0.0000i
-0.4186 - 1.0013i  -0.2851 - 1.7655i
 0.5507 - 0.7589i   1.0478 - 0.4655i

Last 3 terms:

 left channel       right channel
 0.7388 - 0.3032i   0.4740 - 0.6654i
 0.5507 + 0.7589i   1.0478 + 0.4655i
-0.4186 + 1.0013i  -0.2851 + 1.7655i

Looks conjugate symmetrical

I wonder if you are processing the left and right channels at the same time and it is thinking it is a 2d signal or a complex valued signal.

The only way you get an asymmetric fft result is if the signal is complex.


Envelope following using FFT bin power -> Low pass effect on the modulation transfer function

Any set of FFT framings of the data is a form of averaging filter or FIR filter, and thus will have a low-pass filtering effect over some envelope frequency range.  You can understand this effect by treating the FFT window as a FIR filter kernel and calculating its frequency response.

If the width of your FFT window is somewhere near the period of 5 Hz, then you will be in an area of filter roll-off.  You could try a shorter window, compensate for the window roll-off, perhaps using a "flat-top" window if your frequency of interest is well within the window's passband, or use another form of enveloper follower which does not have an intrinsic roll-off in the area of interest.


How do I explain a complex exponential intuitively?

From a graphical point-of view, it is  an infinite spring, whose distance between adjacent coils reflects the frequency of the complex exponential:

If you have a 1D time x-axis, you may be used to draw functions along a single 2nd y-axis  dimension: sines, cosines, etc. If you want to plot a complex function,  you need one x-axis, and 2 y-axes for the real and imaginary parts. You can plot this into a 3D visualization, see above, and you see a spring (although it is some hot summer here). Changing the frequency dilates or contracts the spring. It is known as the Heyser corkscrew, or Heyser spiral, in DSP.

A little more details can be found there:

3D wiggle plot for an analytic signal: Heyser corkscrew/spiral
DSP Concepts Visually Explained

The Fourier transform tells you that any wire (a function) can be reproduced by a superposition of scaled and shifted springs.

After the what, the why. Complex exponentials (or cisoids) are special in  that if one is filtered (with a moving average) it keeps the same shape. So, they are invariant, under  Linear-Time-Invariant (LTI) systems. Invariant vectors/functions are often an appropriate way to study systems or transformations. Morever, as they form an orthogonal basis, they form a basis of choice to decompose arbitrary vectors, and to study how the latter are affected by LTI systems.

Last, the complex exponential is itself invariant under differentiation ($(e^{z})' = e^{z}$), a specific linear and invariant operator), which makes it quite unique, with interesting properties.


Why Does the Odd Multiple of $ \frac{\pi}{4} $ on Gaussian Cause Loss in Repeatability Under Image Rotations?

The answer boils down to 2 issues with the practical approximations of the Gaussian Kernel:

Though the Gaussian Kernel is radially symmetric its discrete approximation has a rectangle support. Unless this support will have infinite length a rotation by any angle different from a multiplication of 90 degrees will yield a shape which has to modified to fit a rectangle. Either choice (Cropping or extending) will yield different results as the weights will be different from the original approximation.
One must pay attention that the longer the support of the filter the less this issue is pronounced.
If I remember correctly, the SURF algorithm uses Box Blur based approximation of the Gaussian Kernel. This approximation, though can be very close visually, doesn't retain the radial symmetry of the Gaussian Blur hence even with large length ha different results when rotated.

I hope this clarifies the issue.


Discrete Fourier Transform for beginners

Yes, m and n are the dimensions of the image.  F(u,v) will be a complex-valued (real and imaginary) function of u and v.
Yes your equation is correct and F(u,v) is the result of the forward DFT. It is correct that F(u,v) has real and imaginary parts, and the multiplication of that into the complex sinusoid will result in a quantity that also has real and imaginary parts.

So how do you get back a real valued image if each of the numbers in the sum over u and v is complex? The key is that F(u,v) has certain symmetry properties due to f(i,k) being real valued.  This symmetry results in the sum over u and v adding up to a real number.

The same thing happens with the one-dimensional DFT.  In that case a real valued signal f(t) produces a conjugate symmetric DFT F(u) (F(u) = F*(N-u)). The symmetry is a bit more complicated in 2D, but that is why this works out.

I often see students manipulate the DFT coefficients and accidentally destroy the conjugate symmetry. Then when they take the inverse DFT the result is complex valued when it should have been real valued.

Note that numerical errors will sometimes produce an inverse DFT that has a very small imaginary part, even though the symmetry is correctly preserved.  In these cases it is best to just take the real part as your result.


In the Fourier domain, the frequency increases with the distance from origo. So an ideal low-pass filter in the Fourier domain is shaped like a disk centered at the origo.

Here they work with only one quadrant of the Fourier domain. Therefore, they create a pie-shaped mask H for the quadrant by computing the Euclidean distance from the quadrant corner and store the distance values in D. The mask H is then chosen as all elements in D with values less than the parameter P (which is the radius of the disk).

So the answer to your question(s) is that the distance computation helps the method compute the filter shape. The key is that a constant Euclidean distance from a point defines the shape of a circle, and all elements inside that circle form the disk that is used as filter mask.

So, naturally, after filtering, only the low-frequency components of the image are left in the Fourier domain. Then the low-pass filtered image is retrieved by inverse transforming the filtered image back from the Fourier domain to the spatial domain.

EDIT

Here is an illustrated example to clarify the answer.

A typical way of illustrating filtering in the Fourier domain is to show the input image.

Then we show the Fourier transform of the image.

Here the origo of the Fourier space is in the center. Since the frequencies increase with the distance from origo, the ideal filter is a disk with its center at the origo. Here we see the disk.

The disk radius is the input parameter P. We use the disk to filter the image by multiplying the filter with the transformed image. We then get the filtered transform.

To get the filtered image we take the inverse transform of the filtered result.

So far so good. As I said, this is usually how low-pass filtering in the Fourier domain is explained.

Now what I was trying to explain when talking about quadrants of the Fourier domain, which got a bit confusing, is that the result from the Fourier transform in the MATLAB code provided does not have origo centered as in the second image in the example above (the image showing the unfiltered Fourier transform).

Instead the result from the Fourier transform F=fft2(double(f)) looks like this (we refer to it as the transform image).

This means that since the filter is a disk with its center in origo, the disk must be placed in the corner of the transform image. Since the Fourier transform is cyclic, the filter should wrap around to cover all four corners of the transform image.

Here is where the distance matrix D comes in. It is used to compute the filter mask above. D is computed by taking the Euclidean distance from each corner (using a neat trick with meshgrid). So Dis the Euclidean distance map looking like this,

where blue values represent low distances and red values represent high.

% Top left corner of D
0.00000   1.00000   2.00000   3.00000   4.00000   5.00000
1.00000   1.41421   2.23607   3.16228   4.12311   5.09902
2.00000   2.23607   2.82843   3.60555   4.47214   5.38516
3.00000   3.16228   3.60555   4.24264   5.00000   5.83095
4.00000   4.12311   4.47214   5.00000   5.65685   6.40312
5.00000   5.09902   5.38516   5.83095   6.40312   7.07107

The filter mask can therefore simply be acquired by selecting the elements in D that have a distance value less than or equal to P.

The rest of the code works as in the example described above. That is, the filter mask is multiplied by the transform image, and the result is inverse transformed to get the low-pass filtered image.

I hope this clarifies what role the distance computation has for creating the filter.


Why not stop with the first level decomposition in wavelet?

The number of levels (frequency range represented) is discretionary to the user usually automatically set as to include a limited number of frequencies which is dependant on the lenght of the series and which is representative of most of the signal.

Applying a wavelet function splits the signal in two coefficients (scaling and detail coefficients) for a number of levels which represent a range of frequencies of this split defined by selecting the number of decomposition levels. Selecting the number of decomposition levels is subjective as per the frequencies you may like to extract.


FASTICA algoritm produce independent components does not equal to the dimensionality of EEG signals

ICA can be seen as a dimension reduction technique similar to PCA.  Although PCA consists of the eigenvectors of the covariance matrix i.e. second order statistics, ICA is not restricted and can contain higher order stat. The goal of ICA is to span the signal space in the most statistical independent vectors. So you might get a considerable dimension reduction.


Second Order Transfer Function : Imaginary component in $R$?

HINT:

If you're looking for the zeros of $H(z)$, try $z_0=1$ or $z_0=-1$. If that works, try to figure out why, and derive it yourself. Then, think about what the magnitude and the phase of those zeros might be.

EDIT (in reaction to your comment): If you have a second-order polynomial with real-valued coefficients, then two complex conjugated zeros are not the only option; you can also have two real-valued zeros.


How to determine initial values in Kalman Filter

but the output signal in the paper is much smoother than my output and

I don't know which paper we're referring to, but it's far from uncommon that you'd smoothen your output for plotting – it would be adequate to specify how you did that, but notice that this sadly doesn't mean everyone adheres to that practice.

  vertical axis is different in my case

that might just be a scaling error, and is not necessarily bad. Again, we don't know which paper is about and whether you're working on exactly the same data, but the fact alone that your line is much thinner and looks like it has more points indicates that they applied some smoothing and/or normalization.

  I think that is because of Kalman filter initialization; initial state and predict.

I don't see indications for that. On the contrary, assuming things converge, I'd say your initial state is fine.


What are typical audio waveform synchronization methods?

Crosscorrelation, find the maximum, and use that to offset one file vs the other.


Why is the Fourier transform so important?

This is quite a broad question and it indeed is quite hard to pinpoint why exactly Fourier transforms are important in signal processing. The simplest, hand waving answer one can provide is that it is an extremely powerful mathematical tool that allows you to view your signals in a different domain, inside which several difficult problems become very simple to analyze.

Its ubiquity in nearly every field of engineering and physical sciences, all for different reasons, makes it all the more harder to narrow down a reason. I hope that looking at some of its properties which led to its widespread adoption along with some practical examples and a dash of history might help one to understand its importance.

History:

To understand the importance of the Fourier transform, it is important to step back a little and appreciate the power of the Fourier series put forth by Joseph Fourier. In a nut-shell, any periodic function $g(x)$ integrable on the domain $\mathcal{D}=[-\pi,\pi]$ can be written as an infinite sum of sines and cosines as

$$g(x)=\sum_{k=-\infty}^{\infty}\tau_k e^{\jmath k x}$$
$$\tau_k=\frac{1}{2\pi}\int_{\mathcal{D}}g(x)e^{-\jmath k x}\ dx$$

where $e^{\imath\theta}=\cos(\theta)+\jmath\sin(\theta)$. This idea that a function could be broken down into its constituent frequencies (i.e., into sines and cosines of all frequencies) was a powerful one and forms the backbone of the Fourier transform.

The Fourier transform:

The Fourier transform can be viewed as an extension of the above Fourier series to non-periodic functions. For completeness and for clarity, I'll define the Fourier transform here. If $x(t)$ is a continuous, integrable signal, then its Fourier transform, $X(f)$ is given by

$$X(f)=\int_{\mathbb{R}}x(t)e^{-\jmath 2\pi f t}\ dt,\quad \forall f\in\mathbb{R}$$

and the inverse transform is given by

$$x(t)=\int_{\mathbb{R}}X(f)e^{\jmath 2\pi f t}\ df,\quad \forall t\in\mathbb{R}$$

Importance in signal processing:

First and foremost, a Fourier transform of a signal tells you what frequencies are present in your signal and in what proportions.

  Example: Have you ever noticed that each of your phone's number buttons
  sounds different when you press during a call and that it sounds the same for every phone model? That's because they're each composed of two different sinusoids which can be used to uniquely identify the button. When you use your phone to punch in combinations to navigate a menu, the way that the other party knows what keys you pressed is by doing a Fourier transform of the input and looking at the frequencies present.

Apart from some very useful elementary properties which make the mathematics involved simple, some of the other reasons why it has such a widespread importance in signal processing are:

The magnitude square of the Fourier transform, $\vert X(f)\vert^2$ instantly tells us how much power the signal $x(t)$ has at a particular frequency $f$.
From Parseval's theorem (more generally Plancherel's theorem), we have
$$\int_\mathbb{R}\vert x(t)\vert^2\ dt = \int_\mathbb{R}\vert X(f)\vert^2\ df$$
which means that the total energy in a signal across all time is equal to the total energy in the transform across all frequencies. Thus, the transform is energy preserving.
Convolutions in the time domain are equivalent to multiplications in the frequency domain, i.e., given two signals $x(t)$ and $y(t)$, then if

$$z(t)=x(t)\star y(t)$$
where $\star$ denotes convolution, then the Fourier transform of $z(t)$ is merely

$$Z(f)=X(f)\cdot Y(f)$$

For discrete signals, with the development of efficient FFT algorithms, almost always, it is faster to implement a convolution operation in the frequency domain than in the time domain.
Similar to the convolution operation, cross-correlations are also easily implemented in the frequency domain as $Z(f)=X(f)^*Y(f)$, where $^*$ denotes complex conjugate.
By being able to split signals into their constituent frequencies, one can easily block out certain frequencies selectively by nullifying their contributions.

  Example: If you're a football (soccer) fan, you might've been
  annoyed at the constant drone of the vuvuzelas that pretty much
  drowned all the commentary during the 2010 world cup in South Africa.
  However, the vuvuzela has a constant pitch of ~235Hz which made it
  easy for broadcasters to implement a notch filter to cut-off the
  offending noise.[1]

A shifted (delayed) signal in the time domain manifests as a phase change in the frequency domain. While this falls under the elementary property category, this is a widely used property in practice, especially in imaging and tomography applications,

  Example: When a wave travels through a heterogenous medium, it
  slows down and speeds up according to changes in the speed of wave
  propagation in the medium. So by observing a change in phase from
  what's expected and what's measured, one can infer the excess time
  delay which in turn tells you how much the wave speed has changed in
  the medium. This is of course, a very simplified layman explanation, but
  forms the basis for tomography.

Derivatives of signals (nth derivatives too) can be easily calculated(see 106) using Fourier transforms.

Digital signal processing (DSP) vs. Analog signal processing (ASP)

The theory of Fourier transforms is applicable irrespective of whether the signal is continuous or discrete, as long as it is "nice" and absolutely integrable. So yes, ASP uses Fourier transforms as long as the signals satisfy this criterion. However, it is perhaps more common to talk about Laplace transforms, which is a generalized Fourier transform, in ASP. The Laplace transform is defined as

$$X(s)=\int_{0}^{\infty}x(t)e^{-st}\ dt,\quad \forall s\in\mathbb{C}$$

The advantage is that one is not necessarily confined to "nice signals" as in the Fourier transform, but the transform is valid only within a certain region of convergence. It is widely used in studying/analyzing/designing LC/RC/LCR circuits, which in turn are used in radios/electric guitars, wah-wah pedals, etc.

This is pretty much all I could think of right now, but do note that no amount of writing/explanation can fully capture the true importance of Fourier transforms in signal processing and in science/engineering


Weak sense stationarity sufficient for multitaper spectrum analysis?

The answers to the questions are

Yes WSS is sufficient.

Quasiperiodic

Semiperiodic


Is there an analogue of coherence between data sets related by a non-linear transformation?

Are you interested in coherence's ability to discriminate correlations at different frequencies, or are you interested in any measure of non-linear correlation? Coherence is basically cross-correlation in the frequency domain.

Mutual information can measure non-linear correlations, but says nothing about the form of the correlation. It's also difficult to estimate accurately without large amounts of data (kernel density estimation is a popular tool to help with that right now, though I think cupolas have begun to gain traction).

I don't know that anybody has done so, but you could alternatively use one signal to predict another using an artificial neural network with one hidden layer. Your correlation metric could then be something like an RMSE between predicted and target signals. It feels like a hack, though. With some work you could try to infer something about the relationship type based on the trained weights of the network. I'm not an expert, but I don't know of a formal framework to actually do that. It's also not an established metric, so you have to also do the work to convince people that it's meaningful (much more difficult if you can't cite sources that did it before you).


Convolution is what happens when a signal passes through an LTI system.
$$y(t)=x(t)\star h(t)=\int_{-\infty}^{\color{red}\infty}x(\tau)h(t-\tau)d\tau$$
The keyword in this question is real-time. The associated keyword in the answer would be weather the LTI system is causal or not.

If the system is causal, then $h(t)=0, \ \forall t<0$, and therefore $h(t-\tau)=0,\ \forall \tau>t$ which means
$$y(t)=\int_{-\infty}^{\color{red}t}x(\tau)h(t-\tau)d\tau$$
It means we only need the inputs only from $-\infty$ up to the present time $t$ which is all available for a real-time implementation.

So the convolution can be implemented in real-time if the system we implement is causal.

To conclude and generalize this idea for two signals, if at least one of the signals are right-sided (which is practically the case), then the real-time implementation is possible.

Regarding your other question, it really depends on the delay performance of your system.


As I've mentioned in a comment, time variant signal is a term which is not commonly used in signal processing, and I also couldn't find a clear definition. After browsing a few texts which use the term (usually written by people not specialized in signal processing), I've come to the conclusion that by time variant signals they actually mean non-stationary signals. The latter being a clearly defined term (check the wikipedia page on stationarity). In short, it means that the signal's properties, such as its mean and variance, change over time. And this is indeed the case with EEG signals, just as it is the case with speech signals.


Assume two discrete-time sequences: $x$ (length $L_\mathrm{x}$) and $y$ (length $L_\mathrm{y}$). Furthermore, $L_\mathrm{x} \geq L_\mathrm{y}$. To calculate the discrete-time (linear) convolution of $x$ and $y$, first pad $y$ with $L_\mathrm{x} - L_\mathrm{y}$ zeros, i.e.
$$
y'(n) =
  \begin{cases}
    y & 0 \leq n < L_\mathrm{y}\\
    0 & L_\mathrm{y} < n \leq L_\mathrm{x} - 1
  \end{cases}
$$
Then calculate the convolution of $y'$ and $x$ as described here, for example.


Isolate the non blurred part of foucsed image

In general, this is an image segmentation problem (http://en.wikipedia.org/wiki/Image_segmentation) into which you would be trying to isolate the focused to the non-focused regions of the image.

Optical lenses are equivalent to low pass filters anyway and the effect of a low pass filter on a signal is to smooth it out by limiting the higher frequency content.

In digital image processing, low frequencies correspond to big objects (general form) and higher frequencies correspond to small objects (details).

Therefore, to isolate the focused region, you could be looking for areas that are at a higher standard deviation relative to the blurred areas which would be at a lower standard deviation because of the smoothing operation of the out-of-focus lens.

To distinguish 'blurred' vs 'non-blurred' regions  you can setup a "window" of dimensions $M \times N$ and slide it over your image.

At each position of the window compute the variance of the "patch" values as if they were a one-dimensional vector of $M \times N$ values. This will result in another "image" (or more generally a two-dimensional array) depicting the variation of the standard deviation at each location on your image. You can then collect all of these values and obtain their histogram. The histogram will exhibit a multimodal distribution (For more information please see http://en.wikipedia.org/wiki/Multimodal_distribution ).

In the case of the image in your post, the histogram is likely to exhibit two modes. One corresponding to relatively high values of standard deviation (focused area) and one corresponding to relatively low values of standard deviation (blurred areas).

You can then position a threshold, for example just between the two modes to establish a criterion between your focused and unfocused regions and then select all "patches" (or MxN window areas) whose standard deviation is higher than your threshold. (For more information please see: http://en.wikipedia.org/wiki/Thresholding_%28image_processing%29)

This will isolate the focused parts of the image for you.

Hope this helps.


Playing samples in a loop, without clicks or pops

The popping comes from discontinuities in the signal or its derivative.  Even if you match the starting sample value with the ending one, you will still have a "corner" in most cases.  Imagine a pure sine wave of one frequency at the start and a pure sine wave of another at the end. Discontinuities tend to produce a lot of high frequency noise.  Your second solution using the low pass filter eliminated that noise (it smoothed out the signal), but at the cost of eliminating the signal too.

The simplest solution to this is just to crossfade the end into the start.  Admittedly, at typical playback rates (22KHz) you are looping about 200ms which gives very little time to fade in/out.  At this point, you may be closer to doing an (overlapped) window function.  Doing an abrupt crossfade would be roughly akin to a Tukey window.